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The standard variational derivation of stellar matter structure in the Wigner-Seitz approximation is generalized 
to the finite temperature situation where a wide distribution of different nuclear species can coexist in the same 
density and proton fraction condition, possibly out of /3 -equilibrium. The same theoretical formalism is shown 
to describe on one side the single-nucleus approximation (SNA), currently used in most core collapse supernova 
simulations, and on the other side the nuclear statistical equilibrium (NSE) approach, routinely employed in 
r- and p-process explosive nucleosynthesis problems. In particular we show that in-medium effects have to be 
accounted for in NSE to have a theoretical consistency between the zero and finite temperature modelling. The 
bulk part of these in-medium effects is analytically calculated in the local density approximation and shown to be 
different from a van der Waals excluded volume term. This unified formalism allows controlling quantitatively 
the deviations from the SNA in the different thermodynamic conditions, as well as having a NSE model which 
is reliable at any arbitrarily low value of the temperature, with potential applications for neutron star cooling and 
accretion problems. We present different illustrative results with several mass models and effective interactions, 
showing the importance of accounting for the nuclear species distribution even at temperatures lower than 1 
MeV. 

PACS numbers: 26.50,+x, 26.60.-c 21.65.Mn, 64.10.+h, 


I. INTRODUCTION 

Since the pioneering work of G.Baym and collaborators in 
the early seventies [l], y], the theoretical formalism to vari- 
ationally calculate the equation of state and composition of 
neutron star crusts with cluster degrees of freedom is well 
settled and has been exploited by many authors in the last 
decades ffl- Because of the crystalline structure of low den¬ 
sity neutron star matter, the problem of neutron star structure 
at a given pressure is indeed reduced to the composition of a 
simple Wigner-Seitz cell, composed of a single nucleus im¬ 
mersed in a gas of neutrons and electrons. The ground state of 
the system is then given by a set of coupled variational equa¬ 
tions for the nucleus atomic and baryonic number (and shape, 
if the more exotic pasta phases are included), the volume of 
the cell, and the free neutron density m. 

An alternative formulation within density functional the¬ 
ory was developed at the same time in another seminal pa¬ 
per for neutron star physics by J.Negele and D.Vautherin |8[|. 
This entirely microscopic approach is in principle more ap¬ 
pealing than a cluster model because it allows accounting for 
the polarisation of the neutron and electron gas. More gener¬ 
ally, a microscopic description avoids the artificial distinction 
between clusters and free neutrons, and naturally accounts 
for the interface interaction between them. For this reason, 
and due to the great improvements of the predictive power of 
mean-field energy density functionals in the last decades [sj- 
III , microscopic Hartree-Fock and Hartree-Fock-Bogoliubov 
methods have been widely employed for the computation of 
the neutron star equation of state lfl3l4l9ll . As a consequence 
of this important collective theoretical effort, present uncer¬ 
tainties on the equation of state of neutron star matter at 
zero temperature are essentially limited to the still imperfect 
knowledge of the density dependence of the symmetry en¬ 


ergy JH, which is itself better and better constrained thanks 
to the recent improvements in ab-initio neutron matter calcu¬ 
lations J2ll - l28tl . 

This standard picture was recently challenged in Ref. j29j] , 
where it is argued that the attractive interaction between the 
lattice nuclei and the surrounding free neutrons might induce 
lattice instabilities similar to the case of the ferroelectric phase 
transition in terrestrial metallic alloys. This very promising 
avenue is not explored in the present paper and we restrict to 
the standard picture where the neutrons-nucleus interaction is 
supposed to only lead to a modification of the surface tension. 

Neutron stars being born hot, a natural extension of the 
above mentioned works concern the consideration of finite 
temperature stellar matter, with applications ranging from 
neutron star cooling, accretion in binary systems and dynam¬ 
ics of supernova matter with associated nucleosynthesis prob¬ 
lems. For these applications matter is typically out of )3- 
equilibrium and therefore needs to be considered in a large in¬ 
terval of baryonic densities pg and proton fractions y p . Finite 
temperature mean-field calculations in the Wigner-Seitz cells 
have been largely employed f30l - l35ll . However, because of the 
computational effort associated to these calculations, micro¬ 
scopic modelling of the finite temperature Wigner-Seitz cells 
is not adapted to the large scale calculations needed for su¬ 
pernova simulations, even if some large scale time-dependent 
Hartree-Fock (TDHF) calculations start to be performed f36i - 
[38lll . For this reason, hybrid models with cluster degrees of 
freedom are more appealing to address the finite temperature 
problem. The extension of the Baym et al. compressible liq¬ 
uid drop models to finite temperature was already proposed 
in the eighties |5|] and allowed the elaboration of the famous 
Lattimer-Swesty (LS) ll39ll and Shen @] supernova equation 
of state models, which are still widely used in present super¬ 
nova simulations. 

The problem of both microscopic and liquid-drop models is 
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that they share the so-called Single-Nucleus-Approximation 
(SNA), that is a unique configuration is assumed for each 
(T,pB,y P ) thermodynamic condition, against the very prin¬ 
ciple of statistical mechanics which stipulates that finite tem¬ 
perature corresponds to a mixing of different microstates. In 
particular, in the LS and Shen EoS, besides free nucleons, 
only one kind of light cluster (a particles) and one kind of 
heavy cluster are assumed to exist. The idea is to account 
in an average way for the properties of the statistical cluster 
distribution. The SNA may not affect very strongly thermo¬ 
dynamic properties of matter in the temperature and density 
domains of interest EH, but it has important consequences 
for dynamical processes dependent on reaction rates of spe¬ 
cific nuclei [42, [43] and for the gas-liquid phase transition. 
Therefore, more modern approaches rely on an extended nu¬ 
clear statistical equilibrium (NSE) concept, where the distri¬ 
bution of clusters over, in principle, all mass numbers is taken 
into account and obtained self-consistently under conditions 
of statistical equilibrium E3E3l- Originally, the NSE was in¬ 
troduced to describe the reaction network taking place at the 
end of the evolution of massive stars in red supergiants 0. 
Being very diluted, nuclei interact weakly and are almost not 
modified by the surrounding medium. These conditions nat¬ 
urally lead to the Saha equations. The NSE in the dense and 
hot matter in the core of supernovae was first applied in the 
EoS of Hillebrandt and Wolff Ell¬ 
in recent NSE implementations EsHUlh the interactions 
between a cluster and the surrounding gas is treated in the 
so-called excluded-volume approach. The clusters and the 
gas of light particles do not overlap in space and the clus¬ 
ters binding energy is kept as in the free limit. It is known, 
however, from virial expansion at low density and quantal ap¬ 
proaches ifsH - EHl . that the cluster properties are modified by 
the coexistence with a_gas. In particular, the recent G.Shen et 
al equation of state [581] includes these in-medium effects for 
light particles within a virial expansion insuring the proper 
model independent low-density limit. Moreover, the excluded 
volume treatment of cluster-nucleon interaction is not compat¬ 
ible with microscopic calculations in the Wigner-Seitz cell, 
where cluster properties are naturally modified by the sur¬ 
rounding gas by the density dependence of the self-consistent 
mean-field and the Pauli-blocking effect of occupied single¬ 
particle states. The conceptual difference between the classi¬ 
cal excluded-volume picture and the quantal picture emerging 
from microscopic calculations was discussed in Ref. (590. It 
leads to two different definitions of clusters in dense media, 
namely configuration-space and energy-space clusters, with 
different particle number and energy functionals. Including 
one or the other of the two definitions in a finite temperature 
NSE partition sum will naturally produce differences in the 
observables, even if the total free energy of the Wigner Seitz 
cell entering in the SNA approaches |39,[4y,[58|l does not de¬ 
pend by construction on the cluster definition 115911 . As a con¬ 
sequence, it is not clear if the NSE models have the correct 
limit towards T = 0, where the SNA approximation becomes 
exact. Recent comparisons 0 indicate that huge differences 
exist among the different models even at very low tempera¬ 
ture, suggesting that the zero temperature limit is not fully 


under control. Such an uncontrolled model dependence might 
be an important hindrance to pin down the EoS dependence of 
supernova dynamics EH EH- 

In this paper we develop an analytical unified theoretical 
formalism to describe on one side the single-nucleus approx¬ 
imation (SNA), and on the other side the nuclear statistical 
equilibrium (NSE) approach. To this aim, we map the en¬ 
ergetics and composition of a microscopic Wigner-Seitz cell 
into a model of the same cell with cluster degrees of freedom. 
If a density and isospin dependent modification of the cluster 
surface energy is included, this cluster model can thus exactly 
span the full spectroscopy (ground state and excited states) of 
the extended Thomas-Fermi (ETF) approximation 0 with 
the only uncertainty given by the employed energy density 
functional. 

A variational minimization of the total free energy of the 
Wigner-Seitz cell with respect to the cell composition leads 
to the standard SNA equilibrium approach, at zero as well as 
finite temperature. A complete finite temperature treatment is 
obtained by calculating the partition sum of a system of inde¬ 
pendent cells, leading to a statistical distribution of cells with 
different compositions. NSE equations naturally emerge from 
this treatment, but energy-space clusters are demonstrated to 
be the correct degrees of freedom in order to get a consistent 
treatment towards the zero temperature limit. We also show 
that a cut-off in the cluster density of states has to be applied 
in order to avoid double counting of scattering states. 

The first part of the paper is devoted to zero temperature. 
Section lHAl defines the degrees of freedom and associated en¬ 
ergy functionals used in this work. Section IU bI gives the vari¬ 
ational equations to be solved at zero temperature to get the 
ground state of stellar matter. The non-standard case where p 
equilibrium is not imposed is also considered. This case is not 
physically realistic, but gives the reference zero temperature 
limit of supernova matter, thus guaranteeing the consistency 
of the finite temperature formalism. To maximize the predic¬ 
tive power of the formalism, an experimental nuclear mass 
table is used in section HTCl to predict the composition of the 
neutron star crust, and results are compared to the rich litera¬ 
ture available on this subject. The equation of state is briefly 
addressed in Section III Dl and to conclude the zero tempera¬ 
ture discussion, the issue of phase transitions is analyzed in 
section lHEl We confirm that the constraint of charge neutral¬ 
ity quenches the first order nuclear matter liquid-gas phase 
transition. A residual very narrow transition region exists at 
densities of the order of po/5 — Po/3, depending on the in¬ 
teraction, which physically corresponds to the emergence of 
pasta phases. 

In the second part of the paper, we switch to finite temper¬ 
ature. Section UlI Al gives the derivation of the coupled vari¬ 
ational equations in the SNA approximation, as well as some 
applications in /3-equilibrium. Sections 1 111 Bl II 11 Cl hui Id the 
partition sum of the model in the canonical and in the grand- 
canonical ensemble, leading to the derivation of the general¬ 
ized NSE equations, which are compared to the SNA approxi¬ 
mation in section llllDI The way in which the phenomenology 
of dilute nuclear matter is modified, in stellar matter, by elec¬ 
trons is discussed in section UlI El Finally section HVl gives a 
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summary and conclusions. 


II. ZERO TEMPERATURE STELLAR MATTER 
A. Energy in the Wigner-Seitz cell 

Let us define a zero temperature thermodynamic condition 
for compact star matter as a given value for the baryon den¬ 
sity ps and proton fraction, y p = p p jP b , where p p is the pro¬ 
ton density. Since by definition there is no interaction among 
Wigner-Seitz cells, the total energy density of the system is 
given by: 


e ws(PB,yp) 


= 1 Ews( 0 


lim 

N ^°° LjLi Vws(i) 


(1) 


where E^ s (i) is the total energy (including rest mass contri¬ 
bution) of the i — th WS cell. Each cell consists of Nws neu¬ 
trons and Zws protons in a volume Vws- We make the standard 
simplifying approximation that the cell consists of a single 
cluster with the possible addition of an homogeneous nuclear 
gas, as well as a homogeneous electron gas. This approxima¬ 
tion is inspired by the numerical results of microscopic cal¬ 
culations 111314 1911 . The polarisation of the nuclear gas by the 
cluster is shown to be small by these works. For this reason, 
this effect is generally accounted in cluster models as an in¬ 
medium modification of the surface tension HiH, see be¬ 
low. Concerning the electron gas, self-consistent calculations 
have shown that, because of the high electron incompressibil¬ 
ity, the homogeneous approximation is excellent for all densi¬ 
ties 10. This Wigner-Seitz picture is however not fully real¬ 
istic since it is also well known that at finite temperature light 
clusters can coexist with the single heavy nucleus |63l|64-]. For 
this reason in the Lattimer-Swesty equation of state a par¬ 
ticles are added to the nucleon gas inside the Wigner-Seitz 
ce 110, but interactions between the a’s and the cluster (or 
the gas) are neglected in that model. This coexistence effect 
of heavy and light clusters will be automatically accounted by 
our formalism, because in the NSE model presented in section 
lIHCl the equilibrium configuration will consist in a mixture of 
different WS cells containing clusters of all species. How¬ 
ever, two-body Coulomb and possibly nuclear effects due to 
multiple clusters inside a same cell are out of the scope of 
the present treatment, and each (light or heavy) cluster will be 
associated to its proper WS cell. 

As we shall explicitly work out, for a given set 
(AwsJws-, V\vs) ( Aws = N w s + Z ws ,l ws = N w $ — Z ws ) equi¬ 
librium imposes a unique mass and composition of the clus¬ 
ter and of the gas. Five variables define this mass and com¬ 
position, namely: the cell volume Vy/g, the gas density and 
composition p g = p ng + p pg , y g = p ng - p pg , where p ng (p pg ) 
is the density of neutrons (protons) in the gas, and the neu¬ 
tron N and proton Z numbers associated to the cluster. The 
total energy in the Wigner-Seitz cell is written as E’/// s = 
Z ws m p c 2 +N W s m „c 2 +E WS with: 


E\vs{A,Z,p g ,y g ,p p ) = E vac + Vws (^hm + zl?) +SE. (2) 


Here, EnM{Pg,yg ) is the energy density of homogeneous 
asymmetric nuclear matter, e'^Pe/) is the total energy den¬ 
sity (including rest mass contribution) of a uniform electron 
gas at the density p e i = y p PB imposed by charge neutral¬ 
ity, E vac (A,Z) is the energy of a cluster with Z protons and 
N = A —Z neutrons in the vacuum, and SE is the in-medium 
modification due to the interaction between the cluster and 
the gas. A part of the in-medium correction is given by the 
Coulomb screening by the electron gas, and by the Pauli- 
blocking effect of high energy cluster single particle states 
due to the gas [56i]. This latter effect can be approximately 
accounted for in the local density approximation by simply 
subtracting from the local energy density the contribution of 
the unbound gas states. This local density approach is cer¬ 
tainly insufficient for the in-medium effect on light clusters 
for which more sophisticated approaches have been proposed 
but we expect it to represent the most impor¬ 
tant correction for medium-heavy nuclei. Indeed residual sur¬ 
face terms in that case appear to have only a perturbative ef¬ 
fect 0. 

As in Ref. [5jl] we introduce the left-over bound part of 
the cluster. 



(3) 

(4) 


that we call ”e-cluster”. po(5) and Po p (8) stand for the total 
and, respectively, proton densities of saturated nuclear matter 
of isospin asymmetry 5 = 1— 2po iP /po- Po(<5) may be calcu¬ 
lated as |0] 


(5) 

where K sat is the symmetric nuclear matter incompressibility, 
and L sym , K sym denote the slope and curvature of the symmetry 
energy at (symmetric) saturation po(0). 

In the above expressions the quantity 5 represents the 
asymmetry in the nuclear bulk. It differs from the global 
asymmetry of the nucleus I/A = 1—2 Z/A because of the 
competing effect of the Coulomb interaction and symmetry 
energy, which act in opposite directions in determining the 
difference between the proton and neutron radii. For a nu¬ 
cleus in the vacuum we take the estimation from the droplet 
model ll65h : 


5 = 5 0 = 


1 I 9 Jsym 1 

40 JT73 


( 6 ) 


In this equation, J sym is the symmetry energy per nucleon at 
the saturation density of symmetric matter, Q is the surface 
stiffness coefficient, and a c is the Coulomb parameter. In the 
presence of an external gas of density p g and asymmetry 5„ = 
(Png — Ppg ) /Pg= yg/Pgi the bulk asymmetry defined by eq.© 
is generalized such as to account for the contribution of the gas 
as I 
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5(Pg,Vg) ~ ( : po(S) ) So(Ze ' Ae) + Po{8) 5g ' (?) 

where So(Z e ,A e ) is the asymmetry value given by eq.© if we 
consider only the bound part of the cluster, A = A e , Z = Z e , 

l = Ie- 

For simplicity, in the following variational derivation of 
the equilibrium equations we shall initially assume 8 = I/A 
in eq.© which implies neglecting isospin inhomogeneities. 
However, we will include this effect in Section |HI Cl 

It is interesting to observe that eq. © can now be conve¬ 
niently written as: 


Ews(A,Z,p g ,y g ,p p ) = E e + Vws (ehm + Ed') +8E sur f. (8) 

with E e standing for the in-medium modified cluster energy 
in the e-cluster representation: 

E e (A , Z, p g , y g , p p ) = E mc - Ehm - F 8E Cou iomb > (9) 

Po 

or, alternatively, 

E ws (A,Z,p g ,y gl p p ) = {Vws- Vo)E hm + VwsE r e < i t 

+ E' aC + SE sur f + SEcoulomb, (10) 

where Vo = A/po(5) is the equivalent cluster volume cor¬ 
responding to 8 isospin asymmetry. In this representation, 
that we call "r-cluster” representation [[59], the in-medium ef¬ 
fects only affect the surface properties of the cluster. The in¬ 
medium bulk term apparent in eq.© is here interpreted as 
an excluded volume. At variance with the classical Van der 
Waals model, this ’’excluded volume” is not a simple limita¬ 
tion of the r-space integral of the partition sum, but it directly 
affects the energetics of the Wigner-Seitz cell. 

Once the dominant bulk and Coulomb in-medium effects 
are accounted for by the definition of the e-cluster representa¬ 
tion eq.©, the residual in-medium energy shift 8E sur f can be 
shown to behave as a surface term [5~9[ roll . 

The different contributions to the energy are defined as fol¬ 
lows. The presence of electrons in the cell modifies the clus¬ 
ter energy with respect to its vacuum value by the so called 
Coulomb shift, 

E»uc = E vac + 5Ecouigmb ^ (n) 

with 

SEcoulomb = a c fwsA~ l ^Z 2 , (12) 

and the Coulomb screening function in the Wigner-Seitz ap¬ 
proximation given by, 

^) ,/ 3 -Kra> u3) 

where we used for the average proton density inside the nu¬ 
cleus pop (5)=po(5)(l —5)/2. For the electron total energy 


density (containing the rest mass contribution) we use the ex¬ 
pression proposed in Ref. 10 ] and valid above 10 4 g-cm 3 
where electrons may be considered free. 


(2f 2 +1) t (r +1) 1/2 - In (r + (r +1) 1/2 )], 

(14) 

where t = h(3n 2 Pel) 1 ' 12 /m e \c. The total electron chemical po¬ 
tential (including the rest mass contribution) is defined as a 
function of the total proton density p p = y p PB as 


m ei c 


/j ptot 

=p,)= 


3 


m„,c 


8(3p f /7T 2 ) 2 / 3 /C 


(15) 


( f 2 + l)V 2 (l +6f 2 ) + f2 ( 2f2 + 1 ) 


(f 2 + l)!/2 (1 T t 2 )l/ 2 _ 


Unless otherwise explicitly mentioned, we will use for the 
energy functional of the cluster in vacuum, E vac (A,Z), the ta¬ 
ble of experimental masses of Audi et al. i?], publicly avail¬ 
able in electronic format. When these latter are not known, 
which is typically the case close and above the drip-lines, we 
will use a liquid-drop parameterization l68ll with coefficients 
fitted out of HF calculations using the same Skyrme effec¬ 
tive interaction which is employed for the homogeneous gas. 
This parameterization, hereafter called Skyrme-LDM model, 
reads: 


I/vac / r )7\ 2 7 2 

- J j^-=a v -a s A~ 1/3 -a a (A) fl - — J -a c -^. (16) 
with the asymmetry energy coefficient: 

a a {A) = -(17) 

1 -L v 


For the numerical applications concerning the NSE model in 
section IIIICI this parameterization will be supplemented in 
the case of even mass nuclei with a simple phenomenologi¬ 
cal pairing term, A (A) = ±12/VA where +(-) corresponds to 
even-even (odd-odd) nuclei. The in-medium surface correc¬ 
tion 8E sur f(A,8,p gl 8 g ) due to the interaction with the sur¬ 
rounding gas can in principle be accounted for by a density 
dependent modification of the surface and symmetry-surface 
coefficients. A determination of these coefficients within the 
extended Thomas-Fermi model [66] will be published else¬ 
where [69]. For the numerical applications of this paper, we 
will ignore this correction, and consider that the main in¬ 
medium effect is given by the bulk nuclear and Coulomb bind¬ 
ing energy shift. 

Below saturation, the Coulomb screening effect of the elec¬ 
trons is never total. This implies that only a finite number 
Nspecies of nuclear species (A,Z) can exist at zero temperature, 
and consequently a finite number of WS cells Nws(PB,yp ) = 
Hspecies- Eq. © then becomes: 


Ews(PB,y P ) = 


N ws(PB,y P ) p (k\ 

V -1 EwsW . . 

(18) 

hm 

N k ,v ^°o V 

(19) 
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where V is the total volume and p(k ) is the number of realiza¬ 
tions of the k-cell. The Single Nucleus Approximation (SNA) 
lil] consists in considering Nws(PB,y P ) = l,p(l) = 1. This 
approximation is exact at zero temperature in the absence of 
phase transitions, and in principle should fail at finite tempera¬ 
ture, even in the absence of phase transitions. In the following 
we explore if phase transitions are there or not, and the degree 
of violation of SNA at finite temperature. 


B. Zero temperature solution in the SNA 


The variational formalism to obtain the composition of 
stellar matter at zero temperature has been proposed long 
ago Cl II and regularly employed since then, using more so¬ 
phisticated models for the nuclear energetics B70lT73[l . We 
will follow the very same strategy, but at variance with the 
seminal papers Cl [2[] , we will determine the optimal con¬ 
figuration for each given (ps,y p ) point without implementing 
-equilibrium in the variational constraints. This choice will 
allow us keeping the same formalism for neutron star crust 
and the finite temperature supernova matter, which is not in 
/3-equilibrium. Concerning the specific application to the NS 
crust, we will determine in a second step the y p (pB) relation 
imposed by /3-equilibrium. 

The variables to be variationally found are 
(VVs,A,5,p ? ,y^). The two constraints, which will lead 
to the introduction of two chemical potentials, can be written 
as: 

Pg = 7 T— (A W s - A e ), (20) 

*WS 

yg = 7 }— (Av.s -4) • (21) 

VWS 

Using the relations ©,(|4} between r-clusters and e-clusters 
we can write the auxiliary function to be minimized: 


@(A,8,p g ,y g ,V W s) = £HM+E e /V W s 

- ap g (po —A/Vws) + «po ( Pb~A/V W s ) 

- PygiPo-A/Vws) + /2poPfl(l -2 y p ) 

A8 

~ PPotj —, ( 22 ) 

Vws 

where a and are Lagrange multipliers. An additional com¬ 
plication comes from the fact that at zero temperature the gas 
can only be a pure gas. 

Indeed within the neutron and proton drip-lines a pure nu¬ 
cleus solution is by definition more bound than a solution 
where one particle would be in the gas. The drip-lines are 
defined by the lowest N (Z) solution of the equations: 

E nuc (N+ l,Z,p e/ ) -E nuc (N,Z,p el ) > 0; 

E" UC (N ,Z + l,p e /) —E nuc (N,Z,p e i) > 0. (23) 


Notice that because of the electron screening the drip-lines in 
the neutron star crust are displaced with respect to nuclei in 
the vacuum, and in particular the fission instability line does 
not exist. However eas. (l23l > admit a solution for any N, Z, 


meaning that when the equilibrium solution is below that line 
the nucleus will be in equilibrium with the vacuum. Above 
the neutron (proton) drip-line, we will have an equilibrium 
with a neutron (proton) vacuum gas. This T = 0 anomaly is 
very well known in nuclear matter. An equilibrium with the 
vacuum does not impose an equality between two chemical 
potentials, because the vacuum has p = 0. If a system with A 
particles and energy E(A) is in equilibrium with the vacuum, 
its chemical potential is defined by a one-sided Maxwell con¬ 
struction between E{A = 0) = 0 and E(A) with slope E (A)/A. 
The chemical potential of this particular equilibrium is given 
by: 



(24) 


which implies d(E/A)/dA = 0, that is the equilibrium solu¬ 
tion minimizes the energy per particle (and not the total en¬ 
ergy, as it is the case in a finite well defined volume V). 

Coming back to the minimization of the auxiliary function 
ea. (l22i l. the minimization with respect to the gas densities 
gives the definition of the neutron (proton) chemical potential 
Pn (Pp) 

a + /3 = — = — if p ng > 0; (25) 

Po Po 

a-1 3 = ^ = if Ppg > 0) (26) 

Po Po 

with p n ( p -) g = d£HM/dPn(p)g- ^ one °f the two densities is 
zero, that is below the corresponding drip-line, we lose one 
equation but also one unknown variable, and one can use one 
of the conservation equations to determine the missing vari¬ 
ables. The result is a system of four coupled equations: 


PBp(n) ~ 


-4(1 T 5) 

2 Vws 


PBn(p) ~ Pg 

g(E nuc /A) n 

d A ls =°’ 

1^1 1 , 

a as U -h 


1- 


PqVws 


A(l±g) 
2 Vws 


1 dE e 


A 38 


1 =F <5 dA 
„ Pg d Po 
%o 2 dd ' 


Is = ±Pg 


1=F<5 


l-Pl 

Po 


(27) 

(28) 

(29) 

+ 

(30) 


The upper (lower) sign refers to a neutron (proton) gas, Pb p („) 
indicates the proton (neutron) baryon density for a neutron 
(proton) gas, p g = p n ( p ) g , Pg = p n { P ) g ■ Th e two last equa¬ 
tions suppose that E nuc is a differentiable function of A and 
8, which is obviously not the case if we take experimental 
masses. In this case the derivatives have to be interpreted as 
finite differences. We can see that the Coulomb screening ef¬ 
fect of the electrons enters in the equilibrium equations, while 
the kinetic energy of the electrons does not play any role in 
the equilibrium sharing. This is the reason why this term is 
usually disregarded out of /j equilibrium. However we will 
see that it does play a role, determining the possible existence 
of phase transitions. 
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Different observations are in order. 

First, from eq.(l29l> we can see that, both below and above 
the drip, the minimization conditions correspond to the min¬ 
imization of the energy per nucleon with respect to the nu¬ 
cleus size, at the isospin value imposed by the constraint and 
the chemical equilibrium with the gas eo. (l30l) . Concerning 
eq.(l30l>. the coupling of the isoscalar to the isovector sector 
is trivially due to the fact that we are using (A,5 = 7/A) as 
isoscalar and isovector variables instead of (A,I) which would 
be the more natural choice if we did not have po = po(<5). 
With this choice of variables, if we consider the textbook ex¬ 
ample of two ideal gases composed of two different species 
of molecules 1, 2, E = E\(A\,I\) +E 2 {Ai 1 h) fulfilling the 
conservation equations 

A = A\ +A 2 ; / = I\ (31) 


using the same Lagrange multiplier method as before and 
defining fi , /I 3 as the conjugated chemical potentials of com¬ 
ponent 2, we find the classical equality of chemical potentials 
if we work with the variables (A,/): 


dEi 

dx; w 


d £ 2 1 dE\ dE2 , 

<9A 2 |/2-A1, dh ,Al “ dh |a2_M3 ’ 


(32) 


while we have a coupling to the isovector sector for the mass 
sharing equation if we work with (A, 8): 


dE 1 
dA\ 


1 x ^E\ 

jU + SijU 3 ; ^r-Ui =Aip 3 - 


(33) 




FIG. 1. Surface of the energy per baryon £jvs/Avrs at p B = 10 ~ 7 
fm 3 and different values of proton fractions. The cluster energy is 
calculated according to FRDM (a) and using the experimental 
database Si(b). 


The solution is particularly simple in the case of symmetric 
nuclei 1 = 0: 


This is a very natural result, because with this choice of vari¬ 
ables the two constraints are not independent any more, that 
is the constraint associated to the /) multiplier contains the 
variable A 1 . 

Second, the factor (1 — p g /po) introduces a coupling be¬ 
tween the two subsystems cluster and gas, that is an interac¬ 
tion. This comes from that fact that our two systems are in fact 
coupled: the energy of the cluster depends on the composition 
of the gas as it can be seen from eq.®. This coupling is due to 
the fact that a part of the high density part of the Wigner-Seitz 
cell is constituted by the gas. From eq.dTOl) , we can see that 
this is an effect of the excluded volume which enters the mass 
conservation constraint. In the e-cluster language (see eq.®), 
we can equivalently say that it is an effect of the self-energy 
shift of the e-cluster inside the gas. This shows that the ex¬ 
cluded volume indeed acts as an interaction. This effect goes 
in the direction of reducing the effective chemical potential 
with respect to the non-interaction case, that is reducing the 
cluster size. If we account for the cluster compressibility, that 
is the 8 dependence of po, an extra effective coupling emerges 
(last term in the r.h.s. of eo. (l30l >). 

In the case of moderate asymmetries below the neutron 
drip, the set of coupled equations eqs. ©.(EE®,® re¬ 
duces to the single equation eq.(l29l> giving the most stable iso¬ 
tope for a given asymmetry. If we assume a functional form as 
eo.(IT6l> for the cluster energy functional, this equation admits 
an analytical solution: 


A eq {I = 0) 


2 a s 

a c { 1 - fws ) 


(35) 


In the vacuum p e i = 0-.fws = 0 and we get a nucleus around 
A ss 55, A eq = 2a s /a c , while at saturation density p e i = po p = 
Po/2 A eq —> °°, showing that we do obtain the homogeneous 
matter limit at saturation. 


C. The structure of the neutron star crust 

The different solutions of eos. (l27l >- (l3()l > lead to a unique 
composition ( A,8,p g ,Vws ) for each couple of external con¬ 
straints (pB,yp). Let us consider the energy density, 
e ws(PB,yp) = Ews/Vws- It contains a baryonic e B and a lep- 
tonic e^ 1 part, e^ s {p B ,yp) = {p p m p + p„m n )c 2 + e B + elf, 
with e B {p Bl y p ) = E e /V W s + Chm (see eq.®). In the long- 
lived neutron star, the proton and neutron densities do not cor¬ 
respond to separate conservation laws because weak processes 
transforming protons into neutrons are in complete equilib¬ 
rium. The structure of the neutron star crust is then obtained 
by choosing, among all the different Z wg values correspond¬ 
ing to different values of the proton fraction y p , the one lead¬ 
ing to an absolute minimum of the energy density. This mini¬ 
mization condition reads: 

£ws(Pb) = min ( £ws(PB,y P )) ■ (36) 

yp 


a a(A) i 2 z 2 

«v-:—— = 2a e (l — Jws)—r- 


a 


A 2 


In the inner crust above the neutron drip the densities are con¬ 
tinuous variables and the energy density is a differentiable 
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function. The minimization then trivially gives the usual 
chemical -equilibrium condition 


P)p tot Pjptot 
OL WS _ ot WS 

dp„ dp p 


2 . de B 2 d £ B 

m n c + - - m p c - -— 

vPn &Pp 


d_e% 

d Pp 


— n tot — n ,ot — u tot — 0 

— Pn Pp Pel ~ U > 


(37) 


where p\ ot = den/ Pp t + m,c 2 with i = ri, p is the chemical 
potential including the rest mass contribution. 

Below the drip-line (outer crust), the baryonic energy den¬ 
sity is simply given by e B = E nuc /Vws =E e {p s = 0) / Vws > and 
the minimization condition (l36t reduces to: 


£ ws(Pb) = mm 


Vws 


-f Pn>n„C -\- PplYlpC~ - 


(38) 


The value of Z leading to the minimal energy Z = Zp eq (A) 
corresponds to the equilibrium nucleus. This is still a /3 - 
equilibrium condition, but it has to be interpreted as the en¬ 
semble of two inequalities: 

Pn\N— 1,Z+ 1) — pf (N, Z) - plf (Z) < 0, (39) 
PT(N,Z)~p t p ot (N+l,Z- 1) - p'Jf (Z- 1) > 0. (40) 

In this set of inequalities, pj"' (N,Z)dp n = + 

1 ,Z) - £ b {N,Z) with <o B (N,Z) = ( E nuc {N,Z ) + Nm n c 2 + 
Zm p c 2 )/Vws(N,Z). An equivalent relation holds for p 1 ^ 1 . 

Most of the published studies on the composition of the 
neutron star crust employ empirical mass formulas A II 
or microscopic functionals from self-consistent HFB calcu¬ 
lations [ 16fc 1 70 h 73II to describe the cluster energy functional 
E vac . A functional approach is unavoidable in the inner crust, 
because no experimental measurement exists above the drip¬ 
line. Conversely, in the outer crust the predictive power of the 
approach entirely depends on the quality of the mass formula 
to describe experimental data. Now, it comes out that the en¬ 
ergy surface in the presence of the electron gas has a huge 
number of quasi-degenerate minima. 

This is shown for an arbitrary chosen density within the 
outer crust p B = 10 7 fm 3 in Figure Q] This figure shows 
the energy surface of the equilibrium Wigner-Seitz cells cor¬ 
responding to different y p values obtained using the FRDM 
parameterization by Moller and Nix from ref. |74| as well as 
the measured nuclear masses from ref. ||67|| . We can see that, 
though the FRDM predictions are very close to the measured 
mass, the obvious tiny differences can affect the determina¬ 
tion of the absolute minimum. This means that even modern 
highly predictive mass formulas describing nuclear masses 
within 0.5 MeV or even less can lead to inexact results when 
applied to the outer crust. 

The importance of this model dependence is shown in Fig¬ 
ure [2] and Table I. 

Let us first discuss the outer crust, on the left of the verti¬ 
cal lines in FigO We can see that the use of an experimental 
mass table leads to sizable differences even with sophisticated 
mass formulas like the FRDM model iil. The solution of 
the variational equations for densities up to about p B = 10 5 
fm -3 , leads to an equilibrium nucleus whose mass has been 




P B ( fm ' 3 ) 


FIG. 2. (Color online). Outer crust composition at T = 0: Baryonic 
(top panel) and atomic (lower panel) numbers of the /3-equilibrium 
nucleus as a function of the baryonic density. BPS corresponds 
to predictions by BPS 0]; FRDM and exp+SLY4(SKMs) stand for 
present model predictions when nuclear masses are calculated ac¬ 
cording to Finite-Range Droplet Model of Ref.j74i] and, respec¬ 
tively, atomic mass data of Ref. t67il + LDM-SLY4 (SKMs) model of 
Ref. 16811 . The vertical lines mark the drip-line in the stellar medium. 


experimentally measured. This means that, up to that den¬ 
sity, fully model independent results can be obtained using 
the experimental mass table, as it is done for the values noted 
as ”exp+Sly4” in the table. If the FRDM model is used in¬ 
stead (results labelled "FRDM” in the table), differences ap¬ 
pear even at very low density. Not only the density at which 
the transition from a nuclear species to another is observed is 
not correctly reproduced (column 1 and 4, respectively), but 
the isotope (column 2 and 5) and element (column 3 and 6) se¬ 
quence is not correct either. These differences are due to the 
imperfect reproduction of nuclear mass measurements by the 
model, and stress the importance of using experimental values 
for the nuclear mass when studying the crust composition. 

In the inner depths of the outer crust, the equilibrium nu¬ 
cleus is so neutron rich that no mass measurement exists. 
Then the crust composition depends on the theoretical cluster 
functional employed, and more specifically on its properties in 
the isovector channel, which are still largely unknown. As it 
is well known, this induces a strong model dependence on the 
composition. As shown in Table I, the lowest density at which 
this model dependence appears is of the order of p B = 10 5 
fm" 3 . At that density the solution of the variational equa¬ 
tions solved using the SLY4 functional when experimental 
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TABLE I. Composition of the outer layer of the outer crust of a cold 
neutron star as a function of baryonic density. FRDM and exp+SLY4 
stand for model predictions when nuclear masses are calculated ac¬ 
cording to Finite-Range Droplet Model of Ref. IH and, respectively, 
atomic mass data of Ref.|i67] + SLY4 model of Ref. {68}] . The atomic 
and mass numbers in italics for the set exp+SLY4 correspond to nu¬ 
clides for which experimental mass evaluations (or extrapolations) 
do not exist. 



FRDM 


exp+SLY4 


Pb (fm 

- J ) 

A 

Z 

p B (fm~ J ) 

A 

Z 

1 .000- 

10 ~ 1U 

56 

26 

1 .000- 10~ lu 

56 

26 

4.467- 

10~ 9 

52 

24 

5.012-10^ 9 

62 

28 

3.388- 

10~ 8 

62 

28 

1.513 -10~ 7 

58 

26 

4.786- 

10~ 8 

58 

26 

1.698-10' 7 

64 

28 

7.585- 

10~ 8 

54 

24 

8.128-10' 7 

66 

28 

1.950- 

10- 7 

64 

28 

9.772-10' 7 

86 

36 

5.623- 

10- 7 

66 

28 

1.862 -10 6 

84 

34 

1.479- 

10~ 6 

82 

34 

6.761 -10- 6 

82 

32 

2.291- 

10~ 6 

84 

34 

1.096-10~ 5 

96 

34 

4.786- 

10~ 6 

80 

32 

6.918 -10 5 

102 

36 

6.761- 

10~ 6 

82 

32 

9.120-10- 5 

104 

36 

1.349- 

10~ 5 

80 

30 

1.148 -10 4 

106 

36 

3.090- 

10- 5 

78 

28 




7.943- 

10- 5 

124 

42 




1.097- 

10~ 4 

122 

40 




1.585- 

10~ 4 

120 

38 





masses are not available, produces as preferred isotope 96 Se 
(A = 96,Z = 34). Now, the smallest Z for which an experi¬ 
mental mass exists for A=96 is Z=35 showing that this solu¬ 
tion is due to a mismatch between the prediction of the SLY4 
functional and the experimental data. The results in italic in 
Table I are therefore not reliable. The observed deviation in 
Fig[2] between exp+SLY4 and exp+SKMs is similarly due to 
the fact that the mismatch is bigger with the less performing 
SKMs functional. The full model independence of the outer 
crust composition is confirmed by the fact thaUrur results for 
the outer crust are in agreement with refs. EDS. This essen¬ 
tially shows that our variational equations are correctly solved. 
In ref. El the model-independent region is slightly larger than 
in our work, because they have used the FRDM model to com¬ 
plement the experimental information when unavailable, and 
this latter, as we have already stressed, has a smaller mismatch 
with experimental data. 

Whatever the predictive power of the mass model, a model 
dependence is unavoidable in the inner crust, where the equa¬ 
tion of state of the pure neutron gas directly enters in the min¬ 
imization equations To illustrate this point, we show 

in Fig{3] the total composition of the neutron star crust ob¬ 
tained with different models. Whatever the equation of state, 
the predictions of eas. (l27l) - (f30l) show that, in the inner crust, 
the mass and charge of the unique nucleus of the WS cell con- 
tinously increase with baryonic density and then suddenly fall 
to zero. The abrupt cluster disappearence occurs because, de¬ 
pending of the employed interaction, at a density of the order 
of po/5 —po/3 homogeneous matter becomes energetically 
more favorable than clusterized matter. 

The precise value of the density corresponding to clus¬ 




P B ( fm ' 3 ) 


FIG. 3. (Color online). Crust composition at T = 0: Baryonic 
(top panel) and atomic (lower panel) numbers of the equilibrium nu¬ 
cleus as a function of the baryonic density. NV stands for predictions 
by Negele and Vautherin 1; BPS+BBP corresponds to predictions 
by BPsOh an d BBPj^l ; exp+SLY4(SKMs) stand for present model 
predictions when nuclear masses are calculated according to atomic 
mass data of Ref. (H + LDM-SLY4 (SKMs) model of Ref.®]. The 
inset in the bottom panel depicts the evolution with baryonic density 
of the total proton fraction. 


ter dissolution depends on the effective interaction mainly 
through the L sym parameter of the equation of state 1181 but 
also on the exact prescription for the cluster surface tension, 
particularly its isospin dependence which cannot be unam- 
bigously extracted from mean field calculations ElEHl- A 
more sophisticated expression for the surface symmetry en¬ 
ergy, different from the one of ref. | 68] was variationally calcu¬ 
lated in ref. Ell for some selected Skyrme models, and slightly 
higher transition densities are consequently reported. 

Because of the abrupt behavior shown by Fig|3] the crust- 
core transition was typically considered as (weakly) first order 
in the literature i] . For this reason the density of cluster melt¬ 
ing is still known in the literature as the ’’transition density”. 
It is however nowadays well established that at the density 
of nuclei dissolution non-spherical pasta can be energetically 
favored, making the transition continuous from the thermo¬ 
dynamic point of view. We will come back to this point in 
section UTEl 

The effect of the nuclear matter equation of state in the pre¬ 
diction of the composition of the inner emst has been studied 
in detail in the recent years It leads to the differ¬ 

ence in Fig0 between the dotted and dash-dotted line, which 
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represent two characteristic equations of state. A more exten¬ 
sive study of the different Skyrme interactions is beyond the 
scope of this paper, however some extra results on this sub¬ 
ject can be found in ref. 1525]. It is interesting to notice that, 
at variance with A and Z, y p (p p ) plotted in the insert of Fig. 
©b) shows no sensitivity to the equation of state. This means 
that the energetics of electrons dominates over the details of 
the nucleon-nucleon interaction. The significant difference in 
atomic number between our model (irrespective of the effec¬ 
tive interaction) and the original Negele-Vautherin Hartree- 
Fock calculation HI is due to the fact that our cluster model 
with the energy density functional [[680 cq. ( I I 6| i contains only 
the smooth part of the cluster energy. The neglected shell ef¬ 
fects are responsible of the emergence of the magic number 
Z = 40 in the Hartree-Fock calculation. We can see that the 
knowledge of shell closures for extremely neutron rich nu¬ 
clei is much more important for the description of the inner 
crust than the isovector equation of state, and it is clear that 
to be predictive, the model at zero temperature should be aug¬ 
mented of realistic proton shell effects, as it is done in the 
Strutinsky approximation by S.Goriely and collaborators 0 
at the obvious price of a greatly increased numerical effort. 
This limitation of ea. (f]~6l > will however not be a serious prob¬ 
lem for the finite temperature applications for which the model 
has been conceived, and which will be studied in the second 
part of this paper. 

The most striking feature of Fig{3] is the huge qualitative 
discrepancy at high density with the original inner crust BBP 
model 0. To understand the origin of this difference Figj4] 
displays the behavior as a function of the baryonic density of 
the mass of the energy-cluster |59i] from eq.©. SKMs |79j] 
and SLY4 0 effective interactions have been considered. We 
can see that the difference between BBP and our approach 
starts when the e-cluster size starts to depart from the r-cluster 
size, that is when the contribution of the neutron gas becomes 
important. In this situation one can expect a modification of 
the surface energy of the cluster according to eq.©. In BBP, 
the in-medium modified surface energy is assumed to be a 
monotonically decreasing function of the gas density, inde¬ 
pendent of the isospin, exactly vanishing when the density of 
the gas reaches the density of the cluster 0. In the language 
of the present paper, this happens when A e = 0 (see eq.©). In 
such a condition BBP clusters are liquid drops with bulk only, 
and their size naturally diverges. However, this approach ne¬ 
glects the energy cost of the isospin jump at the cluster-gas 
interface. It is shown in refs. [52;, [66], in the framework of the 
extended Thomas-Fermi theory, that the in-medium correc¬ 
tion to the surface energy shows a complex dependence on the 
isospin, and specifically behave very differently in symmetric 
nuclear matter and in the equilibrium with a pure neutron gas. 
Only in the case of symmetric nuclei immersed in a symmetric 
gas the transition to the homogeneous core can be seen as the 
simple vanishing of surface energy with diverging size of the 
nuclei; conversely, in the case of (3 -equilibrium matter, the in¬ 
clusion of in-medium surface effects leads to a weak decrease 
of the average cluster size and a slightly advanced dissolution 
of clusters in the dense matter. 

Finally, it is important to stress that results at densities 
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FIG. 4. (Color online). Behavior of the cluster mass number as a 
function of the baryonic density in the inner crust using two different 
Skyrme functionals (SKMs |7 l J| and SLY4 (§]) for both the free neu¬ 
tron energy density and the nuclear masses according to the LDM- 
Skyrme model of Ref. [ 68). The total mass of the cluster is compared 
to the bound part of the cluster, obtained by simply subtracting the 
number of free neutrons according to eq.©. 

higher than about one fifth of normal nuclear matter density 
are not reliable in any of the presented models because of the 
lack of deformation degrees of freedom which could allow the 
appearance of pasta phases 01 - 


D. Equation of state 

A quantity of primary importance when discussing the 
sensitivity of stellar matter energetics to the details of the 
nucleon-nucleon interaction or linking nuclear parameters 
with astronomical observables is the equation of state and, in 
particular, the total energy density - total pressure dependence. 

The total energy density and pressure of the WS cell are 
plotted in Fig. [5] as a function of baryonic density, in com¬ 
parison with the result from the macroscopic BBP 0 and 
the microscopic Negele-Vautherin 0 model. We can see that 
the quantitative value of the energy density obviously depends 
on the model, and more specifically on the effective interac¬ 
tion, but in all cases over the considered density range the 
energy density surface is convex. This means that there is 
no way to minimize the system energy by state mixing, such 
that the system is thermodynamically stable. The discontin¬ 
uous change of the crust composition due to the shell effects 
only leads to very tiny backbendings in the baryonic pressure 
as shown in the insert of Fig© and already observed by dif¬ 
ferent authors 00lzl. These structures can be formally 
interpreted as phase transitions, but are so small that are not 
expected to have any thermodynamic consequence and can 
simply be understood as an interface effect. 

The absence of a phase coexistence region covering a broad 
density domain, well known in the astrophysical context, is 
surprising from the nuclear physics viewpoint because it is in 
clear contrast with the phenomenology of pure baryonic mat¬ 
ter, which is dominated at sub-saturation densities by the nu- 
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FIG. 5. (Color online). Total (baryonic+leptonic) energy density 
and pressure as a function of baryonic density corresponding to the 
neutron star crust (T=0, j6-equilibrium). Experimental J67] andLDM 
binding energies |68] have been used in the outer and, respectively, 
inner crust. The employed nuclear effective interaction are SLY4 | 
and SKMs J79ll . Present results are confronted with those of NV 
and BBP (If 


clear liquid-gas phase transition^0J. One may wonder if this 
difference is due to the fact that we are limiting our analysis 
to a limited part of the two-dimensional baryon density space 
that is explored in /3-equilibrium. Indeed the j3 -equilibrium 
trajectory corresponds to very neutron rich matter, and it is 
well known that the coexistence zone in the nuclear matter 
phase diagram shrinks with increasing asymmetry. We there¬ 
fore turn to demonstrate that the difference between stellar 
matter and nuclear matter thermodynamics is not restrained to 
/3-equilibrium. 


E. Phase transitions in the inner crust? 

In the previous section we have assumed that a one-to-one 
correspondence exists between baryonic density and chemical 
potential, that is a unique Wigner-Seitz configuration can be 
systematically associated to each pressure and chemical po¬ 
tential field inside the star. 

This is only correct in the absence of phase transitions, and 
it is in principle possible that a mixture of different Wigner- 
Seitz configurations might lead to a lower energy density than 
a periodic repetition of the same cell. The highly degener¬ 
ate energy minima showed by the experimental energy sur¬ 
face even without the (more model dependent) inclusion of 


unbound neutrons beyond the drip-lines (see figure Q]) evoke 
the possibility that first order phase transitions could even ap¬ 
pear at finite temperature in the outer crust. 

More generally, it is well known that such a first-order 
phase transition covers almost the whole phase diagram of 
sub-saturation neutral nuclear matter [80] and has baryonic 
density as order parameter. It is therefore natural to ask 
whether such a phase transition persists in the stellar context. 
As a matter of fact, the existence of such first order phase 
transition is systematically assumed in most seminal papers 
on the stellar matter equation of state [2, ;5j], and in particular 
it is implemented in the publicly available and popularly used 
LS tables f39f] ■ Even more modern equations of state of su¬ 
pernova matteij 50j, jUt] invoke the persistence of the nuclear 
liquid-gas phase transition in the stellar context, based on the 
fact that the baryonic energy density of star matter is unstable 
with respect both to thermodynamic [8(1 HUl and to finite size 
density fluctuations 0i£3, [84]. 

On the other side, it was shown in different works that the 
liquid gas phase transition in stellar matter is quenched by the 
very strong incompressibility of the electron background m 
1821 [84148611 . and microscopic modelling of the Wigner-Seitz 
cell has confirmed a continuous transition from the solid crust 
to the liquid core through a sequence of inhomogeneous pasta 
phases El00 ,[H(tO,[H4H(H- 

It is therefore important to examine this question in further 
detail. 

We have already seen in section III Cl that the solution of 
eas. (l27l >- (l30t is always unique, even if many different solu¬ 
tions can be very close in energy per nucleon (see FigureQ}. 

This means that at zero temperature a unique cluster-gas 
configuration can be associated to a given value of p^ s = 
Aws/Vws , y^ S = Zws/Aws- This statement can of course 
be model dependent, as we have seen that very small varia¬ 
tions of the mass functional can lead to very different results. 
However, even if multiple solutions of the cluster configu¬ 
ration would occur (which will indeed be the case at finite 
temperature), this would not lead to a first order phase tran¬ 
sition. For a first-order phase transition to occur, solutions 
corresponding to different densities should be degenerate in 
(constrained) energy. Then, the absolute energy density mini¬ 
mum would be obtained by mixing these degenerate configu¬ 
rations with different p^ s ,y^ s - If this was the physical result, 
the single-nucleus approximation would fail, and even at zero 
temperature one should account for a distribution of different 
Wigner-Seitz cells. 

Thermodynamic instabilities and eventual phase transitions 
in systems with more than one component have in principle 
to be studied in the full N-dimensional density space [89|]. In 
our case this means that the energy density has to be studied in 
the full two dimensional (p,,.p p ) plane, and the /I-equilibrium 
condition has to be applied only after the Gibbs construction 
is performed (indeed p -equilibrium has to be imposed only 
at the macroscopic level, and can very well be violated at 
the microscopic level of a single cell). However the prob¬ 
lem simplifies if the order parameter is known. In that case it 
is useful to introduce a Legendre transformation of the ther¬ 
modynamic potential with respect to the chemical potentials 
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of all the densities except the order parameter iflcfl . Then 
the multi-dimensional Gibbs construction exactly reduces to 
a one-dimensional Maxwell construction on the residual den¬ 
sity. 

In the case of stellar matter the neutrality condition p p = p e i 
allows a variable change (p„,Pp) (Pb = Pn +Pp,pL = Pel)- 
Due to the very huge electron incompressibility it is reason¬ 
able to expect that the two coexisting phases, if any, would not 
present any jump in electron density [84j]. Microscopic calcu¬ 
lations [* 14i] have convincingly shown that the electron polar¬ 
ization by the proton distribution is negligible, as long as the 
clusters have linear dimensions of the order of the femtome- 
ter. Then, we can safely perform a Legendre transform with 
respect to p/ and introduce the constrained energy density 

£ws(Pb,Pl) = £ws~ PlPl, (41) 

where Ews{pB,pL ) = E ws /V W s, Pl stands for lepton chemi¬ 
cal potential and pi is the value taken by the lepton density 
at chemical potential Lip, p/ = pi(pi). Note that LL/ = 0 
corresponds to /3-equilibrium in lack of neutrinos, p' 0? = 
p' p l + p'Jf . We do not therefore need to examine the whole 
Pl plane, but can limit ourselves to the single point Pl = 0. 

We can then conclude that we can identify the possible pres¬ 
ence of phase transitions in the neutron star crust by simply 
considering the Pb density behavior of the energy density in 
/3-equilibrium, £ws(Pb, Pl(Pb, Pl = 0)). As in section HTbI 
Eb = £ws ~ is obtained solving, for each condition ( PB,yp ) 
the coupled equations (I27l> - (f30l >. 

To better evidence possible convexities, it is useful to intro¬ 
duce a linear bias with slope A#: 

£\vsp B (Pb,Pl = 0) = £ws(Pb ,Pl = 0) — A bPb ■ (42) 

For each A b value, which plays the role of an external chem¬ 
ical potential field, the equilibrium density of star matter is 
given by the minimum of this function. If the function Ews 
is convex, it will be characterized by a single mimimum value 
giving the usual relation between intensive and extensive vari¬ 
ables 

(43) 

oPb 

In this equation, we have introduced a prime symbol on the 
chemical potential to indicate that the electron contribution is 
included in Ews- However, if Ews has concave region(s) on the 
baryonic density axis ps, it will be possible to find one or more 
values of A b such that two (or more) different configurations 
correspond to the same value of the constrained energy den¬ 
sity. This will indicate a first order phase transition, and the 
associated A b value will correspond to the transition chemical 
potential. The constrained energy density eq. (l42l) for cluster- 
ized matter in the crust is displayed in Fig. [6]for some chosen 
values of A b corresponding to minima in the outer (a) and in¬ 
ner (b), respectively. We can see that both in the outer and in¬ 
ner crust the constrained energy surface is smooth and that the 
equilibrium configuration is given by a single Wigner-Seitz 
cell, thus justifying our variational procedure. The FRDM 


-2 

x 10 




FIG. 6. (Color online). Constrained energy densities as a function 
of density, (a) Outer crust as obtained using FRDM. The numbers 
accompanying the curves are in MeV and stand for A#, (b) Crust 
(thick line) and j8-equilibrated homogeneous matter (thin line) at T=0 
corresponding to SLY4 and A^ = 11.25 MeV. 


mass model is limited to nuclei below the dripline and can¬ 
not be used for calculations in the inner crust. For this reason 
we have switched to the Skyrme-LDM mass model [68] to 
produce the right panel. Again, a unique clusterized solution 
characterizes the equilibrium up to a chemical potential of the 
order of 10 MeV (A b = 11.25 MeV for Sly4). At that point, 
the corresponding equilibrium density is of the order of po/3. 
As we have already discussed commenting Figj3j the precise 
value depends on the EoS and on the surface tension. We 
can see from Fig. [6]that at this transition value of the chemi¬ 
cal potential the constrained energy density of the clusterized 
system is equal to the one of homogeneous matter, meaning 
that it is possible to put in equilibrium the two phases. 

This defines a tiny region of first-order phase transition, 
much less extended than the liquid-gas phase transition of 
normal nuclear matter. Indeed this latter covers the whole 
sub-saturation density region. Moreover, we believe that this 
residual phase transition might be an artifact of the present 
model which does not account for deformation degrees of 
freedom. It is well known that in this density domain de¬ 
formed pasta structures have to be accounted for @]. For this 
reason, we do not perform any Gibbs construction and simply 
put to zero the cluster mass at the transition point, assuming 
that pasta would take over. The results of Figj6]show that the 
thermodynamics of /3 — equilibrated matter is completely dif¬ 
ferent from the one of nuclear matter. 

As previously discussed in refs. f84] within mean-field ar- 
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FIG. 7. (Color online). Constrained energy density (eq. ( 144b ) in the 
thermodynamic condition corresponding to the liquid-gas phase tran¬ 
sition of symmetric matter given by Xg = —15.97 MeV and A 3 = 0. 
The different curves represent different values of the total proton den¬ 
sity p p . Solid (dashed) curves correspond to constrained baryonic 
(baryonic+leptonic) energy density. The considered effective inter¬ 
action is SLY4. 


guments, this difference is due to the huge electron gas in¬ 
compressibility which quenches the phase transition in stellar 
matter. 

To demonstrate this point in the framework of the present 
model, we turn to consider the behavior of the baryonic part 
of the energy density £3 = £\vs — £ ‘ e f in the full (p/i-P,,) plane. 

To better spot convexities in the two-dimensional space, we 
introduce again a constrained energy density 

hsMiPS’Pp) = £ h(,Pb-P p) ~ hiPB - h(Pn ~ 2p p ); (44) 

where A# and A 3 represent an isoscalar and isovector external 
chemical potential field. 

Again, phase transitions will be signalled by the existence 
of one or more values (A#, A3) such that two (or more) differ¬ 
ent configurations correspond to the same value of the con¬ 
strained energy density. 

In the case of uncharged nuclear matter, we know that the 
dominant part of the plane is characterized by con¬ 

cavities. It is therefore not surprising to see that this is clearly 
the case for the £g function plotted in Fig. [^corresponding to 
A B = -15.97 MeV and A 3 = 0. 

This figure displays the energy density obtained solving the 
variational equations eas. (127H3()b . biased by an external chem¬ 
ical potential field according to cq. (l44l i. The value of the total 
proton density p„ is the same in each point of the different 
curves plotted in the figure. Because of the neutrality con¬ 
straint p p = p e [, each curve represents a given screening fac¬ 
tor to the cluster Coulomb energy according to eci.( f]~3l . The 
minimum of each curve then gives the ensemble of optimal 
Wigner-Seitz cells corresponding to the chosen A3 value, and 
to different baryon density. The absolute minimum corre¬ 
sponds to the equilibrium Wigner-Seitz cell associated to the 
couple (As, A3). The [N.Z) sequence of the corresponding 
cluster gives the composition, as a function of baryonic den¬ 
sity, of matter at that A3 chemical potential. 


The choice A 3 = 0 selects the equilibrium solutions for sym¬ 
metric matter. A unique point is the absolute minimum for all 
choices of A b except A b = —15.97 MeV which is shown in the 
figure. At that chemical potential, if the electron part of the en¬ 
ergy density is not taken into account (solid curves), two dif¬ 
ferent points correspond to the same constrained energy. This 
corresponds to the well-known nuclear matter phase transi¬ 
tion which at zero temperature takes place at a chemical po¬ 
tential equal to the saturation energy, Pb = —15.97 MeV for 
the SLY4 functional chosen in Fig. [3] We can notice that the 
low density phase, which is predicted to be the vacuum phase 
in mean-field calculations, is obtained here at a low but finite 
density, corresponding to the most stable N = Z isotope 56 Ni. 
This is due to the limitation of mean-field calculations that do 
not account for clusterization at low density. 

In the stellar matter case, because of the Coulomb coupling 
between protons and electrons, the lepton part of the energy 
density is not independent of the baryon part. This means that 
the energy in the Wigner-Seitz cell has to include the electron 
zero point energy as written in eq.© . The total energy den¬ 
sities are given by dashed curves in Fig0 This contribution 
is a simple constant shift of each curve because of the con¬ 
dition p p = p e i, and therefore does not change the sequence 
of optimal compositions as a function of the density. From a 
thermodynamic point of view, we can say |! 86 f] that the canon¬ 
ical solution is the same as without the electron contribution. 
However, the electron energy density is a monotonically in¬ 
creasing function of p e i = p p , and the optimal p p monotoni¬ 
cally increases with pn in this symmetric matter situation we 
are considering. As a consequence, no value of A b can be 
found such that two different Wigner-Seitz cells can be put in 
equilibrium, and the phase transition disappears. This can be 
easily understood mathematically considering that the optimal 
energy density gains an extra term as 

£b -4 £ws = £b + elf (p P ) ■ (45) 

The relations (1431) between density and chemical potential are 
shifted because of the electron contribution 

—>• p' B = Pb+ -p[f ; A 3 —> P 3 = P 3 — -plf , (46) 

and the curvature of the constrained energy density becomes: 

d 2 i ws dp B 1 dp‘ e f 
dpi dp B 2 dp e i 

Because of the very high electron incompressibility, the con¬ 
vexity observed in the baryonic part of the energy density is 
not present any more in the total thermodynamic potential. 
This is known in the literature as the quenching of the phase 
transition due to Coulomb frustration |84ll85]l, and shows | 86 l 
that convexities in the (free) energy density do not necessarily 
correspond to instabilities in the physical system. 

This shows that if one wants to formulate the equilibrium 
problem in the grandcanonical ensemble, one has to account 
for the electron zero point motion. This is a triviality for the 
zero-temperature problem, since the Wigner Seitz cell is natu¬ 
rally defined in the canonical ensemble. However, in the finite 
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temperature NSE problem, which is typically treated grand- 
canonically, this kinetic contribution is usually disregarded 
with the argument that, the electrons being an ideal gas, the 
corresponding partition sum is factorized 115 (X D 10 . It is then 
important to stress that a negative eigenvalue in the baryonic 
energy curvature matrix should not be taken as a sign of a first 
order phase transition: stellar matter inside the convex region 
is clusterized but perfectly stable, and no Maxwell or Gibbs 
constructions should be performed to get the equation of state. 

We will come back on this point in the second part of this 
paper. 

To summarize the results of this section, first order phase 
transitions could be possible in zero temperature stellar mat¬ 
ter only if the function e w $ = Ews/Vws, with E\ys defined in 
eo. lflOt and the values of A,I, p s ,Vws obtained from eas. (l27l> - 
(EoK presents a convexity anomaly as a function of pg, for 
fixed values of p p . As it can be seen in Figs[6] [7] except the 
narrow density domain where deformation degrees of freedom 
have to be accounted for, this is not the case even using energy 
functionals for the nuclear masses which include shell effects. 
In particular, the minimum of the constrained energy density 
for a given set of chemical potential is systematically associ¬ 
ated to a single Wigner-Seitz cell, characterized by a unique 
composition in terms of cluster and gas mass and composi¬ 
tion. This means that the SNA is perfectly adequate to deal 
with the zero temperature problem. 


III. FINITE TEMPERATURE STELLAR MATTER 


In the second part of this paper we extend to the finite tem¬ 
perature regime the modelling of the Wigner-Seitz cell with 
cluster degrees of freedom, presented in the first part. We 
will start by deriving the classical equations corresponding to 
the single nucleus approximation. This approach is at the ori¬ 
gin of most extensively used equations of state for supernova 
mattcr|39. 40. 581. Then the main part of the paper is devoted 
to the derivation of an extended nuclear statistical equilibrium 
model, which by construction reproduces the results of SNA 
if only the most probable cluster is considered, and the same 
chemical potentials are considered. Since the SNA naturally 
converges at T=0 to the standard modelling of the neutron 
star crust, the consistency between the theoretical treatment 
of neutron star crust and finite temperature supernova matter 
will thus be guaranteed. 


A. Single nucleus approximation 


configuration k = ,A^ k \ S^ k \ p g k \y^ } 

F W s{A,Z,p g ,y g ,pp) = F e p - TV WS \nz% M - TV ws ln Z f 

+ SF sur f, (48) 


where (3 = T~ l is the inverse temperature, Zp M (p g ,y g ) is the 
mean-field partition sum for homogeneous matter and Fp rep¬ 
resents the cluster free energy in the Wigner Seitz cell (defined 
below). 

The electron contribution is independent of the different 
configurations and the associated partition sum z e p {p P ) will be 
factorized out. Similar to the previous section, we will neglect 
the surface in-medium corrections to the free energy, though 
they might turn out to be important in the situations where 
the gas contribution is not negligible. We note that in this 
section, at variance with the notations used in Section II, ther¬ 
modynamical quantities corresponding to homogeneous mat¬ 
ter component and electron gas bare the HM and el labels as 
superscripts. 

If we consider temperatures higher than the solid-gas phase 
transition temperature, the free energy of a cluster defined 
by the variable couple (A,Z) or equivalently (A, 8 ) is differ¬ 
ent from its ground state energy because at finite temperature 
the cluster can be found in different translational and internal 
states. 

To calculate this term, one has to consider that within the 
Aws total number of particles, a number .A, = p g Vws belongs 
to the gas part. The entropy associated to these particles is 
already contained in the term In z^ M ■ To avoid a double¬ 
counting of the number of states, the canonical partition sum 
of the cluster must thus be defined summing up the statistical 
weight of the different energy states associated to this reduced 
particle number (see eqs.Q.Q): 


z ws( A ’ s >Pg,y g ) = ex P (~P F p) 


=EE ex p 

p e * 



(49) 



+ E e + E 



The cluster center-of-mass motion is a plane wave. The first 
sum is thus given by the plane wave density of states, with pe¬ 
riodic boundary conditions at the cell borders. This is simply 


E 


Vws 
(27r/i) 3 



(50) 


Notice that the available volume for the center of mass is the 
whole Wigner-Seitz volume, and there is no excluded volume 
effect. The center-of-mass momentum integral is a Gaussian 
integral 


The natural extension at finite temperature of the model pre¬ 
sented in the first part of this paper consists in keeping the 
SNA approach, and replace the variational problem of the en¬ 
ergy density minimization with the variational problem of the 
free energy density minimization. In the e-cluster representa¬ 
tion eq.© the energy is additive and we can write for a given 



The sum over the cluster excited states has to be cut at the 
average particle separation energy, to avoid double counting 
with the gas states. This leads to a temperature dependent 
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degeneracy factor defined by: 


r<S> 

E ex P(~P E *) = / dE*p AS (E*)exp(-pE*) 
e* Jo 

= gp(A,8,p g ,y g ), (52) 

where < S >= min(< S n >,< S p >) is the average particle 
separation energy. For the numerical applications of this pa¬ 
per, we will use for simplicity a different higher energy cut¬ 
off for each cluster species < S >~< S > (A,8,p e /) = min(< 
S„ > (A,Z,p e i),< S p > (A,Z,p e i)), with separation energies 
calculated from the smooth part of the cluster energy func¬ 
tional, given by E]j‘ )M + SEcoulomb■ In the zero temperature 
limit gp —> gcs = 2Jos + 1 gives the spin degeneracy of the 
cluster ground state. 

For the numerical applications of this section, in order to 
be able to study the whole subsaturation density domain in 
equilibrium without any mismatch in the cluster energy func¬ 
tional, we have s yste matically used the Skyrme-LDM model 
for nuclear mass 16811 and considered gos = 1 ■ The level den¬ 
sity p A §(E) is here taken for simplicity with a simple Fermi 
gas formula m. A more realistic choice will be presented in 
section ITlI Cl 

The cluster free energy in the Wigner Seitz cell than reads, 

F$ = Fp* + T]nz$ M (p g ,y g )^ (53) 

= E e (A,8,p g: y g ,p p )-T\nVws-Tlncp - ^T\nA e , 

with Fp uc = Fp ac + 8 F cou i (the equivalent of eq. (fTTV). cp = 

gp{mT /(27r/r)) 3 / 2 , and m the nucleon mass. 

The auxiliary function to be minimized is the extension of 
cq. (f22t including the entropy terms: 


%( A , 8 ,p g ,y g ,Vws) 


FwsiA, 8,p gl y g ,p p ) 

Vws 

a Pg (Po —A/Vws) + ap 0 ( Pb —A/Vws ) 
Pyg (Po ~ A/V ws ) + PqPb(1 -2 y p ) 
PpoA8/V ws . 


The variational equations result: 

dE e | _ Po~ Pg , 

PA ^ B 


P3 


Po 8 -y g 

Po 


3 T ppVws 

2 A poVws — A 
dE e 

pg \ A >pg>yg ~~ P 3 ^' 


Po 
d In cp 

1 PA | 5 -p*-y*> 
d8 po Po Po 


PgVws 


A'* pdp® 

2 dS {po~ p g )(poV W s — A) 


d\ncp 

dS 


W/A) 

l<5,V%S ® 


dA 



d In z% M 

Pb 

= -T 

P 

d Pg 



d\nzP M 

P3 

= T 

d yg 


(55) 


(56) 




(57) 

(58) 

(59) 
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FIG. 8. (Color online). Mass (top panel) and atomic (bottom panel) 
numbers of the unique nucleus of the WS cell as a function of bary- 
onic density for Y p = 0.2 and different temperatures T = 0,1,2.5,5 
MeV. Predictions of present SNA model (solid lines) are compared 
with the LS results 13911 corresponding to LS220 as calculated in Ref. 
j90i] (dotted lines) as well as with the NSE prediction for the most 
probable cluster (dashed lines). The LDM-SKMs model of ref. f68il 
is used for the cluster energy functional. 


The finite temperature predictions of SNA are plotted in 
Figs. m and |9] along predictions of zero-temperature SNA of 
Section II, results of Lattimer-Swesty model lf39ll as available 
in Ref. ^0] and NSE results (see Section III.C). Different 
density, temperature and proton fraction conditions are con¬ 
sidered. The considered effective interactions is SKMs & 
Fig® corresponds to the case where the proton fraction is 
kept constant and equal to 0.2. It shows, as expected, a mono¬ 
tonic decrease of the cluster size as a function of temperature. 
More interesting, the results converge for T —> 0 to our zero 
temperature results in the Wigner-Seitz cell, which we know 
to be exact at the thermodynamic limit, and model indepen¬ 
dent below neutron drip. Fig. [9] illustrates cluster mass and 
charge numbers as a function of baryonic density for different 
temperatures at j3-equilibrium. The observed non-monotonic 
behavior is due to the strong decrease of proton fraction. In¬ 
deed, with increasing density, the proton fraction becomes 
so low that loosely bound hydrogen and helium resonances 
dominate over heavy clusters which dissolve into homoge¬ 
neous matter. At constant proton fraction this effect is not 
apparent, meaning that the in-medium bulk energy shift is not 
enough to suppress the cluster binding. In that case, the pref¬ 
erential cluster size monotonically increases with density up 
to the point where homogeneous matter is energetically pref- 
ered. As in the previous section, we have indicated that point 
by putting to zero the A(pb) and Z(pn ') curves. We cannot 
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FIG. 9. (Color online). The same as in Fig. [ 8 ] for -equilibrium 
and different temperatures T = 0,0.5,2.5 MeV. Predictions of the 
pres ent SNA model (solid lines) are compared with the LS results 
|3fll corresponding to LS220 as calculated in Ref. Hi (dotted lines) 
as well as with the NSE predictions for the most probable cluster 
(dashed lines). The LDM-SKMs model of ref. [ 68 ] is used for the 
cluster energy functional. 


exclude that the inclusion of surface in-medium effects, ne¬ 
glected in this paper, could change this behavior. However 
preliminary results 152fl indicate that this effect is small. 

The qualitative behavior of the cluster size and charge with 
density is similar to the one of the LS220 Lattimer-Swesty 
equation of state, plotted with dotted lines in Figs[ 8 ] and [9] 
Quantitative differences exist nevertheless. On one hand they 
could be due to the (slightly) different equation of state pa¬ 
rameters and cluster surface tension model. Effects of em¬ 
ployed effective interaction on clusters have been already seen 
in Fig. [4]where SKMs and SLY4 have been considered. Prob¬ 
ably more important, the SNA model of Lattimer-Swesty ad¬ 
ditionally accounts for alpha-particles that can be present in 
the WS cell together with heavier clusters. The presence of 
an isospin-symmetric bound component in the gas obviously 
modifies the cluster size and composition. Finally, to obtain 
the emergence (at low density) and dissolution (at high den¬ 
sity) of clusters, first order phase transitions to an a particle 
gas and to homogeneous matter respectively, are implemented 
in the LS mode 1 13911 . We also note that at the highest temper¬ 
atures, our SNA clusters tend to be smaller than in LS. This is 
probably due to the high energy cut in the density of state in¬ 
tegral ea.(l52l> implemented in order to avoid double counting 
of the continuum states, which reduces the statistical weight 
of heavy clusters. We will discuss in section HIlCI that the in¬ 
clusion of the proper statistical weight of clusters of all sizes 


FIG. 10. (Color online). Contours of cluster constrained free energy 
Fp ~ ^bA — A 3 (A — 2Z) (in MeV). The considered thermodynamical 

conditions are: pg = 10 ~ 3 fm~ 3 , Y p =0. 39 and T=0.5 MeV (a) and, 
respectively, 2.4 MeV (b). The values of the external isovector and 
isoscalar chemical potentials (A, A 3 ) are (-10.06 MeV, 3.98 MeV) 
and, respectively, (-11.03 MeV, 5.86 MeV). Experimental values I67il 
and predictions of the 10-parameter mass model of Duflo-Zuker 19 111 
have been used for the binding energies. 

naturally leads to the emergence of an important fraction of 
light particles, and to the disappearance of heavy nuclei in the 
dense medium, without invoking any phase transition. 

In the treatment of finite temperature we have presented in 
this section we have assumed that, similar to the zero temper¬ 
ature case, stellar matter in a given thermodynamic condition 
(PB,yp,T) is characterized by a single well defined Wigner- 
Seitz cell. This is of course an approximation, since what has 
to be minimized at equilibrium is the total free energy density, 
and not the single-cell free energy density. The two varia¬ 
tional approaches will give approximately the same result if 
the free-energy landscape has a single deep minimum. If, 
by contray, different Wigner-Seitz cells correspond to com¬ 
parable free energies, they will all be present in the equilib¬ 
rium configuration, even if with different probabilities. As we 
could have anticipated from inspection of Fig[T] the free en¬ 
ergy landscape is highly degenerate for stellar matter. This 
is confirmed by FigJTO] where, for two arbitrary representa¬ 
tive thermodynamic conditions, ps = 10 3 tin 3 , Y p = 0.39, 
T =0.5 MeV and T =2.4 MeV, the constrained cluster free ener¬ 
gies — A/;,4 — A 3 (A — 2Z) with conveniently choosen values 
for (A b ,A 3 ) are depicted. In each of these plots the different 
nuclei are immersed in the same neutron, proton and elec¬ 
tron gas. At the lowest considered temperature we can see 
that, though a single minimum exists, which corresponds to 
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the solution of the SNA variational coupled equations (l55l >- 
( l59t for this given set of constraints, different (heavy as well 
as light) nuclei might lead to comparable free energy densi¬ 
ties, and will therefore be present at equilibrium. The higher 
considered temperature shows a different pattern. Indeed, a 
plethora of nuclei with masses and isospin symmetries span¬ 
ning important ranges have close values of the constrained 
free energies. It is easy to anticipate that if pairing and shell 
effects were ignored, the minimum of the contrained free en¬ 
ergy would have been even flatter. 


B. Thermodynamic limit in the canonical ensemble 

The very principle of statistical mechanics tells us that 
at non-zero temperature different realizations of the Wigner- 
Seitz cell will be possible within the same constraint of total 
density and proton fraction. 

If we consider a very large volume V which contains a num¬ 
ber n —> °° of different Wigner-Seitz cells for a total number 
of particles A tol and a total isospin asymmetry I tot , a possible 
realization of the system is now characterized by k = {nf j , i = 

1 ,... , 00 } where nf 1 is the number of realizations, within the 
volume V, of an arbitrary Wigner-Seitz cell constituted of a 
cluster with particle numbers * ,/j^ = A < ( !' 1 — 2 and a gas 
with particle numbers Ag’^ = V^ s Pg k \lg’ k ^ = V^\yf\ No¬ 
tice that the gas density and isospin can in principle depend 
on the realization (k) but do not depend on the cell (i). Indeed 
the nucleon gas density is uniform over the volume because 
we have divided it in cells only for convenience, and the vari¬ 
ation of gas particle numbers is just due to the variation of 
Wigner-Seitz volumes. 

The total number of particles in the cell A^ S ,I^ S varies 
from one cell to the other, but the total number of particles 
in the volume V is the same for each realization (k): 


v-E.fW+vWpf), 

i 

(60) 

i 

(61) 

sT 

w 

11 

> 

(62) 


Since we are at the thermodynamic limit, these three condi¬ 
tions are in reality only two 


by the neutrality condition which has to be fulfilled in each 
cell 




V, 


ws 


= P P = y P PB. 


(65) 


We can further simplify the problem considering that, for a 
given macroscopic set of constraints (r,p£,y p ) we will have 
a unique partitioning in the macroscopic volume between the 
cluster fraction and the gas fraction, that is the one that min¬ 
imizes the total free energy. It is very easy to improve on 
this approximation if necessary, by considering the canoni¬ 
cal probability associated to each partitioning. We do not do 
it because it comes out that there are very few combinations 
of p c i = Y,i n iAe ^/V and p ? which lead to the same p P - This 
means that we consider that p g and y g do not depend on (k) 
but only on the macroscopic constraints. Then the conserva¬ 
tion law simplifies to: 


A tot 1 t-i (k) . ( i) , , /*z£\ 

pB = — = -2 rj n) ) Ay+Pg = p cl +Pg, (66) 
hot 1 ^ W A0 1 1 utn\ 

yB = — = ~)_ i n)'r e ’ +y g =y d +yg\ (67) 


and the different realizations of the set of constraints 
( T,p B ,y p ,p g ,y g ) are defined by k = {nf\i =A^\l^} . 

The probability pp (k) of this realization is determined by 
the usual maximization of the information entropy under the 
constraint of the average energy and a sharp constraint on the 
mass A tot and isospin hot eas. d66b .d67i. We define the total Q 
free energy of each realization (k) as: 

Ftot (k) = F cl (k) - TV In (z™ (k)zf) , (68) 

with 

^/W = E»f ) ^(0; (69) 

i 

Fp(i) =E e -T\nV -T\ncp- ^T\nA e . (70) 

It is interesting to remark that the cluster free energy at the 
thermodynamic limit ea. d70l) differs from the cluster free en¬ 
ergy in the Wigner-Seitz cell ea. d53l) . Indeed the number of 
states for the center-of-mass motion has to be calculated over 
the whole volume: 


Pb = 


E ; n| k) (A®+vg> s pP) 


v 

, ho, L-», w (/i° + v^yf) 

yB v v\„Wi/( 0 


E inT% 


ws 


(63) 

(64) 


We can then characterize a realization ( k ) by the fragment 
distribution and the gas isoscalar and isovector densities k = 
{nf 1 ,i= 1 ,..., 00 , Pg k \yf ■*} where now rif 1 is the number of 
occurrences of the Wigner-Seitz cell ( i ) constituted of a gas 
pj, k \yf \ a cluster AIg 11 and a volume V^ s uniquely defined 


E ex p 

p 


a P 2 ] 

V 

f2nmA e \ 

r 2mA e 

(2zr/z) 3 

l J8 ) 


(71) 


This is well known from solid state physics and leads to the 
Bloch theorem: even if the ions are localized at fixed positions 
in the Coulomb lattice, their center of mass motion is a plane 
wave over the whole volume ll92ll93ll . 


1 note that here “total” has another meaning than in Eq. m 
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Thanks to the thermodynamic limit, the partition sums are and 
now factorized 

Zp(PB,y P ) = Y,exp[-PF c i(k)] Up M z e p) . (72) 

k v ' 


In z c a(nA/V,nZ/V) = ——In 
H Vws 


(o e p (A,Z) 


(79) 


The probability of realization (k) is then simply given by: 

Pp (k) = exp [-/3 F d (k)\, (73) 

p 

with 

z p(PB,y P ,P g ,y g ) = £exp[-/3F c/ (k)]. (74) 

k 

Zp is the canonical partition sum of an ensemble of fully 
independent clusters, for a total mass number A c / = V p d and 
isospin I d = Vy d . We note by passing that we can easily ex¬ 
tend this result to the case where we allow mixing of different 
p gl y g by considering cci.(l72b as a constrained partition sum. 
In that case, the total partition sum has to be defined as 

ZpipB-Jp) = Y.Zpi.PB-.yp^P'g-.y'g) (z/j) • 

(75) 

The explicit calculation of Z ^ is a classical problem ll5ll 
[94l] , and its solution is given by: 



in 

(i)A,Z 



where 


(Op(A,Z) = exp 


-f5F e p(A,8,p gl y g ) 


(76) 


(77) 


the sum runs over all possible realizations of the system such 

(k) 

that the total number of particles is A d , n A z is the number of 
occurrences of cluster A,Z in the realization k, and the prod¬ 
uct runs over r-clusters A,Z or e-cluster A e ,Z e equivalently, 
since the two are scaled by a factor which is constant if p g 
and y g are constant. This partition sum can be calculated with 
a Monte-Carlo technique|51] or also analytically via a recur¬ 
sive relation ll86l l94ll . 

Notice that for a finite system the total volume V to t = 

(k) 

Y,i n ) Vi is a fluctuating quantity, and only A d is the same 
event by event. However this is a not problem, because the 
conservation law is applied to the total density. 

It is instructive to consider the SNA limit of a representa¬ 
tive cluster. Let us suppose that the average multiplicity den¬ 
sity <n AZ >_/V k, 1/ < V ws {A,Z,p B ,y B ,p g ,y g ) > for a given 
A = A, Z = Z and < n^z >~ 0, VA ^ A, Z ^ Z, or equivalently 
let us suppose that we consider only the most probable clus¬ 
ter in the partition sum. Since A ,./ = nA, I d = n(A 2Z) we 
immediately get 


Zjj (nA, nZ) 



(78) 


where we have used the Stirling approximation neglecting the 
—n term: ln(«!) ps nlnn — n ps n\nn, and we have introduced 
the free energy densities as — Tln^ = — T\nZp/V. Using 
ea.(T72t and ea.(l48l> the partition sum becomes 

— Tlnzp {pB,y P ) = y Fws(A,Z,p g ,y g ,p p ). (80) 


We can see that we recover a SNA expression which we have 
already shown converges towards the exact result at zero tem¬ 
perature. 

The value of A, Z can be deduced from the equations of state 


„ T dlnz p\ . 

Pb T \y B , 

(81) 

„ T dlnz P 1 

M3_ dy B 

(82) 

which can also be written as 


n d(-Tlnzp-p B p B ) l 

3 p B lyB ’ 

(83) 

d{-T\nzp- piys). 
ay, w 

(84) 

Integrating these equations leads to: 


d (-T\nzp - PbPb ~ k(y B )) =0; 

(85) 

d (-Tlnzp - p 3 y B - h(p B )) =0, 

(86) 

or also 


d(-TInzp -p B p B -piyB) =0. 

(87) 


This is exactly the same minimization problem as in section 
IIII Al with the difference that now the variables are A, 5, Vws 
because p g ,y g are fixed. This physically means that the fact 
of considering a large number of Wigner Seitz cell has elim¬ 
inated the conservation constraint between A, 8 and p g ,y g - 
density and isospin fluctuations are allowed in each Wigner- 
Seitz cell because the conservation law applies only to the 
macroscopic system. 

As a consequence, the equilibrium sharing equations are 
slightly modified: 
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( 88 ) 


U,pj,y s ;(89) 
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with 


lem. Indeed we have: 


Pb = 


(90) 

dpg ’ 

P3 = 

d\nz™ 

d y g 

(91) 


We can also notice that in the limit T —> 0 the sharing equa¬ 
tions at T = 0 that we have obtained, by imposing exact con¬ 
servation laws within the cells, are recovered as they should. 
Indeed in this limit the system is periodic and the global con¬ 
servation law is equivalent to a local (within the cell) conser¬ 
vation law. 


C. the grandcanonical NSE 

The canonical treatment of the previous section is formally 
correct, but has the disadvantage of being extremely expen¬ 
sive from the computational point of view. For this reason, a 
grandcanonical formulation appears more appealing and has 
been preferentially invoked in the star matter literaturej45i- 

Esnamsi. 

To formulate this problem, we consider a very large vol¬ 
ume V —> °° which contains a number n — > °° of a-priori 
different Wigner-Seitz cells, and introduce two external La¬ 
grange multipliers to impose the average isoscalar and isovec¬ 
tor densities over the whole volume. As in the previous sec¬ 
tion, a possible realization of the system is noted by an index 
k = {«. , i = 1 ,..., °o| where rij is the number of occurrences, 
within the volume V, of an arbitrary Wigner-Seitz cell consti¬ 
tuted of a cluster with particle numbers AW,/W and a gas with 
particle numbers A g ^ = p g k \l g ^ = . The total 

number of particles in the cell Ay} s ,I$ s varies from one cell 
to the other, but the total number of particles in the volume V 
(or more precisely, the total density and proton fraction, since 
we are at the thermodynamic limit V —> °°), are fixed by the 
externally imposed chemical potentials p and p 3 . These den¬ 
sities, as well as the average cluster multiplicities < 
and the gas densities < p g >, < y g >, is what we want to cal¬ 
culate. 

As we have already discussed, the gas density and isospin 
could in principle depend on the realization (k) but do not 
depend on the cell (i). The Wigner-Seitz volume is uniquely 
defined by the neutrality condition eq.(l65t in the cell 


pp=E n 


... z^ 

(k) ^e j- Ppg v ws 


V 


(94) 


showing that our variational variables {nf 1 } are not indepen¬ 
dent variables. As it is well known in the framework of the 
self-consistent mean-field theory (95], an equivalent one-body 
problem can be formulated corresponding to the same infor¬ 
mation entropy, therefore to the same set of occupations as 
in the self-consistent problem, but with a different free en¬ 
ergy corresponding to independent particles, which contains 
rearrangement terms. These rearrangement terms will explic¬ 
itly appear in the one-body occupations of the self-consistent 
problem. The free energy of the equivalent one-body problem 
is given by: 


F%{k) = -TVlntfipp) +£»J k) F^(i), (95) 

i 


with 


*a *(0 = 


d F ln 


dn 




= Fp(j) ~ TV^f In (k) + p el (z« + p&Vg 5 ) , 


(96) 


where the derivation is taken at constant p g ,y g ,rij ,j ^ i. The 

(k) 

grand-canonical occupations n] are determined by the free 
energy of the equivalent one-body problem, meaning that they 
directly depend on the electron chemical potential. Notice 
that in principle also (i) depends on the total proton den¬ 
sity through the Coulomb screening term, therefore it should 
also give rise to extra rearrangement terms. However this extra 
term dF^/dp p ■ dp p /dni <x y , is negligible in the thermo¬ 
dynamic limit. The total Gibbs one-body free energy of each 
realization (k) is obtained by Legendre transformation with 
respect to the total baryon number and isospin. This amounts 
to introducing as usual two chemical potentials p^p’^ accord¬ 
ing to: 


G%(k) = F%(k) - Enf ] (p’bA^s + p^ s ) (97) 

i 


We can see that we can define auxiliary chemical potentials as 
Pb = Pb~ Pel/2 ; pi= P 3 + Pei/2 (98) 


w(kk) _ 

V WS ~ 
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(k) ' 
Pp-Ppg 


(92) 


The total Helmholtz free energy of each realization (k) is 
given by eq.(| 


rt«(k) = E »i k) ( F p(i) - TVwP ln ( 4 M ( k ) 4 )) ■ (93 > 


We can see that, because of the dependence on p p of the elec¬ 
tron free energy, this equation defines a self-consistency prob- 


such as to make formally disappear the electron contribution 
in the cluster free energy: 


-16 


(*) = -TVlnzf+E^ (Fp® - TV^hnz H p M (k)) 

i 

_r „WL /HO i v(‘ 4 )1 ,, (A') 1 

/ , n i yPlt A-V ws Pg J + P3 yle + r ws y g J 


(99) 


Using the mean-held relations of uniform nuclear matter: 

ln #2U = ln 4 M ( p «’ y g) + PPBP g + PP3y g - (100) 
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we can see that the gas densities are uniquely determined by 
the external chemical potentials, and independent of the real¬ 
ization, as we could expect: 


Pg = T 
y g = T 


Jl n JIM 
6/111 ~/W 3 

Jl n nHM 
^ 11 U /W3 

<9^3 


IM3 ’ 


I MB' 


( 101 ) 

( 102 ) 


We can then write V^f 1 = vf/ g , = p g , y g k> = y g , and: 


GZ(k) = -TVlnzp 


+E 4 


(*) 


j /3MbM 3 


(4°,/i°) - TV, 


(‘,k) ]n JIM 
WS llu PMflM3 


(103) 


where we have defined the in-medium modified cluster Gibbs 
energy: 


Gpusn 3 (A,S,Pg,y g ,Pp) = Fp — p B A e — p 3 I e . (104) 

The thermodynamic limit implies that all the realizations cor¬ 
respond to the same (infinite) volume: 

V = E (105) 

i 


meaning that the gas contribution becomes completely inde¬ 
pendent of the cluster contribution, and fully determined by 
the chemical potentials 

G) b ot (k) = -TV In 

(106) 

We are ready to calculate the one-body equivalent grand- 
canonical partition sum 

Z 16 = £exp[-/3G£(*)] = (107) 


with 
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(109) 

( 110 ) 


With eg. i ll 1 0b have recovered a NSE-like expression for the 
cluster multiplicities 


< n i Mb,M3 — VpUBHiW 


(HD 


= exp 


~P (Fp(A,8,p g ,y g ,p p )- p B A e -piI t 


( 112 ) 


where the electron energy density and entropy density are 
known. 


It is interesting to notice that the baryonic component (clus¬ 
ters as well as gas) only depend on the baryonic part of the to¬ 
tal chemical potentials p B ,P 3 - These chemical potentials are 
not the thermodynamic potentials conjugated to the densities 
which determine the thermodynamics; indeed they are 
shifted of the electron contribution. This explains why the 
phase transition is quenched in stellar matter even if the bary¬ 
onic chemical potential ll b has a backbending behavior as a 
function of the baryonic density. This point will be further 
discussed in section IIII El The backbending behavior of p B 
was observed in refs. 1811 1861. but it was interpreted as a sign 
of ensemble inequivalence]!^] °r of instability {8 lj], since the 
rearrangement terms coming from the electron contribution 
were not discussed in those papers. 

From a practical point of view, the numerical implementa¬ 
tion of the NSE model is simpler than the one if its approxima¬ 
tion, namely the SNA. Indeed the variational character of the 
approach is fully exhausted by the construction of the partition 
sum, and no extra variation of the energy functional has to be 
performed. This means that we can easily use fully realistic 
functionals for the cluster free-energies with no extra numeri¬ 
cal cost. For the applications shown in this section, we use the 
tables of experimental masses of Audi et al. [67]| and, in or¬ 
der to extend the pool of nuclei for which pairing and shell 
effects are accounted for, evaluated masses of Duflo-Zuker 
JH for the vacuum energies, the full list of low-lying reso¬ 
nances for light nuclei for the degeneracy factor gp, and real¬ 
istic level densities fitted from experimental data from ref. |9(1 
in ea. (l52l >. Only when this information is not available (or 
to make the quantitative comparisons with SNA as in section 
IHIDb . we switch to the Skyrme-LDM mass model. Moreover, 
we explicitly consider isospin inhomogeneities in the spatial 
distribution of clusters due to Coulomb and skin effects. This 
is done considering that the bulk asymmetry 8 entering in the 
in-medium correction to the cluster energies does not coin¬ 
cide with the global asymmetry I/A = 1 — 2Z/A, as proposed 
in ref. [65] (see eq.©). This equation is consistently solved 
with eq.Q which gives the isospin dependence of the satura¬ 
tion density {59]. It was shown in ref. 1 59ll66tl that accounting 
for the difference between bulk and global asymmetry is a 
crucial point to obtain, within a cluster model, energy func¬ 
tionals compatible with microscopic calculations. Again, the 
difference between 8 and I/A is neglected in the numerical 
applications of section llllDI in order to compare the NSE and 
SNA approaches within the same definitions for the physical 
ingredients. 

Our final result, eas. d 10 lt . d 1021 . (1112b is formally very 
close to the different existing versions of grandcanonical ex¬ 
tended NSES® ® ® M !H. This is not surprising, 
since these equations simply state that all the different bary¬ 
onic species are quasi-ideal gases of independent particles. 
However, some specificity of the proposed approach should 
be stressed. 

It is clear from the microscopic treatments of the Wigner 
Seitz cell at zero temperature that any realistic finite tem¬ 
perature model has to include in some way interactions be¬ 
tween the clusters. The way of implementing these in medium 
effects is not unique, and the different treatments lead to 
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a considerable spread in the predictions of extended NSE 
models l60ll . 

The viewpoint we have taken in this paper is that the very 
definition of a Wigner-Seitz cell implies that WS cells are the 
correct variables that can be treated as independent degrees 
of freedom. This fully fixes the in-medium effect under the 
unique hypothesis that each cell contains only one bound clus¬ 
ter. As we have discussed in the introduction, this hypothesis, 
which is employed by all the existing models in the litera¬ 
ture, is certainly not completely correct in general and some 
cluster-cluster interaction should be taken into account! 64] to 
improve the present description. 

The result of building a model on independent WS cells is 
that a NSE-like expression can be recovered for the cluster 
abundances, but with some specific features which insure that 
the zero temperature limit is properly obtained. Specifically, 
we can see from eg. (11 1 2b that the variable conjugated to the 
chemical potentials is not the physical cluster size (A,Z) but 
the reduced value ( A e ,Z e ) (cqs.Q.ijTj) which represents its 
bound part. Moreover, the cluster free energy has to be mod¬ 
ified according to eq.(l53t if one wants that in-medium effects 
are limited to a modification of the surface tension. 


D. NSE versus SNA 


To compare in greater detail the SNA to the NSE results, 
we can evaluate the most probable cluster mass and isospin 
A, I predicted by the NSE. This is obtained by maximizing the 
argument of the exponential in eq. dl 1 111 : 

dGf PM = d ( F p - ^A e ~ Me) =0. (113) 

Since p g ,y g are fixed, we can equivalently put to zero the par¬ 
tial derivatives with respect to A e , I e , or with respect to A, 8 . 
The first choice leads to: 


r)F e 
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dJ 1 A e,Pg,y g ’ 

(114) 

?F e 

1 

dA e |/e,p **' 

(115) 


These equations look very different from the equilibrium 
equations (f55l) . (f56t corresponding to the SNA. However they 
are far from being in a closed form. Indeed, the dependence 
on A e ,I e of G e is highly non trivial: 


Fp(A,8) = Fp(A(A e , 8(A e ,I e )),8(A e ,I e )), (116) 


where the dependence of 8 on A e ,I e is obtained from the so¬ 
lution of the two coupled equations 


A e = A 
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y s 
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(117) 

(118) 


This coupling will induce an effective coupling between the 
isoscalar and isovector chemical potential. After some algebra 



FIG. 11. (Color online). Constrained cluster free energy for p# = 
10~ 3 fm~ 3 , T = 1.5 MeV in j8-equilibrium with chemical potentials 
corresponding to the NSE (full line) and SNA (dashed line) model, 
for a fixed cluster proton fraction corresponding to the minimum of 
the constrained free energy. The arrow gives the SNA solution. The 
dotted line gives the NSE multiplicity distribution in arbitrary units. 
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( 120 ) 


These equations are similar, but not identical to the SNA equa¬ 
tions (I55l) . (f56l >. The difference arises from the fact that the 
Wigner-Seitz volume as a variational variable in the SNA 
approach induces a complex coupling between the different 
equations. In the NSE, the most probable Wigner-Seitz vol¬ 
ume is trivially defined by the condition 


Vws 


Z e 

Pp ~ Ppg 


( 121 ) 


Conversely, we have seen that the same result as in this 
grandcanonical approach is obtained if we consider a canoni¬ 
cal problem with a large number of Wigner-Seitz cells. This is 
not surprising, because the neighboring cells act as a particle 
bath. 

This result implies that we do not necessarily expect that the 
most probable cluster obtained in the complete NSE model 
exactly coincides with the result of the SNA approximation. 
This effect however turns out to be very small. A much more 
important source of difference between SNA and NSE is ex¬ 
pected when the NSE distribution has multiple peaks of com¬ 
parable height. In that case the ( PB,yp ) of the total distribution 
for a given set of chemical potentials is not the same as the one 
of the most probable cluster. This induces a non negligible 
shift between SNA and NSE even at very low temperatures. 

This point is explained in Fig(TT] which shows the in¬ 
medium modified cluster free energy ea. (l53l > as a function of 
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the cluster size. The free energy has been constrained with 
two Lagrange multipliers, corresponding to the chemical po¬ 
tentials obtained in the SNA and NSE model at an arbitrarily 
chosen thermodynamic point belonging to the j3 -equilibrium 
trajectory, p# = 10~ 3 fm~ 3 , T = 1.5 MeV, y p = 0.08503. To 
allow a one-dimensional representation, a cut has been done 
with a plane whose (1—2 • Z/A)-value is as close as possible 
to the corresponding value of the constrained free energy min¬ 
imum. The observed staggering stems from discrete values 
of (1 — 2 • Z/A) and, for A <16, structure effects accounted 
for in the experimental binding energy. The NSE abundan¬ 
ces are also represented in arbitrary units. We can see that 
the NSE abundancies correctly follow the shape of the con¬ 
strained free energy, as implied by eq. dl 121) . This means that 
for identical values of the chemical potentials in the two mod¬ 
els, the optimal SNA cluster (indicated by an arrow in the ex¬ 
ample shown in the figure) should exactly coincide with the 
most probable NSE cluster. However, allowing clusters of 
any arbitrary size and composition obviously alters the map¬ 
ping between density and chemical potential. The deviation 
of chemical potentials is typically very small (for the example 
shown in the figure, we have p# = —12.929(—12.889) MeV , 
p 3 = 13.547(13.507) MeV for NSE (SNA)), but it is sufficient 
to modify the position of the constrained energy minimum. As 
a consequence, a SNA treatment cannot correctly identify the 
most probable cluster. 

NSE predictions were already compared to the SNA 
approximation, both at fixed proton fraction and in /3- 
equilibrium, in Figs. [8]and[9]above. We have seen that, except 
the very low densities where light cluster degrees of freedom 
are important, at low temperature the NSE model is very close 
to SNA. However the consideration of clusters of all sizes nat¬ 
urally leads to a reduction of the cluster size at high density 
and high temperature, similar to the LS equation of state be¬ 
cause of the particular treatment of a particles in that model. 
It is however important to notice that in the complete NSE a 
particles are only abundant for matter close to isospin sym¬ 
metry, while more neutron rich hydrogen and helium isotopes 
prevail in neutron rich matter. This aspect, which by construc¬ 
tion cannot be addressed in the LS model, will be discussed in 
greater detail later. 

At higher temperature the NSE distribution is spread over 
a large domain of cluster sizes and isospin (see the panel (b) 
of Fig. m and the deviation both with SNA and with the 
LS equation of state becomes very large. In particular, the 
abundances are dominated by light resonances and the heavy 
cluster yield becomes increasingly negligible with increasing 
temperature. 

A more detailed comparison between SNA and NSE is 
given by Figs. [12] and [13] in terms of the unique/most prob¬ 
able cluster mass and, respectively, relative mass fraction of 
unbound nucleons. For NSE, the most abundant cluster mass 
is plotted against the average mass of heavy clusters arbitrarily 
defined as clusters with A > 20. 

As in the previous figure. Fig. [12] shows that, whatever the 
density, increasing temperature leads to an increased deviation 
between the average SNA composition and the most probable 
NSE cluster. A huge part of this difference can be explained 



T (MeV) T (MeV) 

FIG. 12. (Color online). Structure of the (P)NS crust at 
/3—equilibrium as a function of temperature for different values of 
the baryonic density ps = 10 6 , 1CL 4 ,10~ 3 ,1CL 2 firD 3 . Solid and 
dashed lines correspond to predictions of SNA and, respectively, 
most abundant cluster in NSE. Dotted lines correspond to the average 
heavy (A > 20) cluster atomic mass in NSE. Standard NSE predic¬ 
tions are plotted agaist predictions of a modified NSE (NSEm) where 
no cluster lighter than A = 20 is allowed to exist. In SNA Skyrme- 
LDM [681 binding energies have been used. In NSE, experimental 
data [67] have been used for the binding energies of nuclides with 
A < 16 and Skyrme-LDM (68tl predictions otherwise. The consid¬ 
ered effective interaction is SLY4. 


by the importance of accounting for (a variety of) light clus¬ 
ters which are entropically favored at increasing temperature. 
This is confirmed by red curves that correspond to a ’’mod¬ 
ified” NSE obtained by artificially switching to zero the sta¬ 
tistical weight of all clusters lighter than A = 20. We can see 
that neglecting light clusters considerably approaches SNA to 
NSE, even if residual differences still persist partially because 
of the shift in chemical potentials discussed above. A com¬ 
plementary view is offered Fig. [13] At intermediate densities 
(p B = 10 4 ,10~ 3 fm- 3 ) and T > 0.5 MeV, SNA and NSE 
predictions agree in the percentage of unbound nucleons, thus 
indicating that the chemical potentials of the two models have 
close values (recall that at a given temperature the gas density 
only depends on the chemical potential). 

At variance with this, the extreme densities show a strong 
reduction of the nucleon gas in NSE with increasing tem¬ 
perature. At the lowest density displayed = 10 6 fm~ 3 
where matter as a whole is close to isospin symmetry, this 
comes from the enhanced production of 2 H and 4 He at the 
cost of unbound nucleons. At the highest density Pb = 10 2 
fm -3 the opposite holds. The extreme neutron enrichment 
of /j-equilibrated matter favors copious production of isospin 
asymmetric hydrogen and helium isotopes leaving thus less 
unbound nucleons. The suppresion of clusters with A < 20 
(red curves) confirms the above reasoning by showing a per¬ 
fect agreement between SNA and NSE everywhere except 
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FIG. 13. (Color online). The same as in figure[l2]for the unbound 
nucleons mass fraction. 


Pb = 10 6 fm 3 and T >1.3 MeV where p c i/pB — > 0 . 

The global behavior of the /) -equilibrated matter composi¬ 
tion in the NSE model is shown in Figs. |T4l[T5] In Fig. fl4lav- 
erage mass, charge and mass fraction of heavy (A > 20) nuclei 
are plotted as a function of density for temperatures ranging 
from 0.4 to 2 MeV. Fig. [15] presents, for the same temper¬ 
atures, the mass fractions of unbound neutrons and protons 
together with the mass fraction of different light species ( 2 H, 
J H, 4 He, / '- 4 H, - 5 He) as a function of density. As men¬ 
tioned in the figure captions, in these cases experimental val¬ 
ues Ei and DZ10 [91] predictions have been used for nuclear 
masses. The unbound nucleon component is treated according 
to SLY4. 

We can again observe the nice convergence towards the 
zero temperature composition of the Wigner-Seitz cell, as well 
as the complex behavior as a function of density for all tem¬ 
peratures, leading to a melting of the clusters in the nuclear 
medium at a density of the order of pp = 0.01 fm~ 3 . As it 
is well known in the literature, the exact value of the transi¬ 
tion density depends on the effective interaction. We do not 
try to make such a study here because the presence of defor¬ 
mation degrees of freedom in the form of pasta phases, here 
neglected, would most probably modify the value of the tran¬ 
sition density. Inspection of Fig. fl5l reveals the importance 
of accounting for all the different light nuclear species and 
not limiting to deuteron and a-particles. This is true for any 
proton fraction, but particularly clear in the very neutron rich 
matter implied by /)-equilibrium, where light unbound reso¬ 
nances completely dominate, together with unbound neutrons, 
the matter composition at high temperature. The considera¬ 
tion of light particles of all species, including heavy hydro¬ 
gens and helions, is natural and easy in the context of a NSE 
model. However to our knowledge no SNA approach includes 
such particles in the description of the average Wigner-Seitz 
cell, even if a very promising step in this direction was re¬ 
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FIG. 14. (Color online). NSE results at )3—equilibrium for differ¬ 
ent densities and temperatures (expressed in MeV and listed in the 
key legend). Upper and middle panels: average mass and atomic 
numbers of clusters heavier than A = 20. Lower panel: correspond¬ 
ing mass fraction. Experimental and DZ10 j9lll nuclear masses have 
been considered for clusters. The SLY4 effective interaction was 
used for the unbound nucleon component. 


cently undertaken in refs.|63 ! , 64]. This underlines again the 
importance of going beyond the SNA approximation in the 
finite temperature stellar problem. 


E. Thermodynamics and electrons 

The problem of the grandcanonical formulation which has 
been recently observed f8ll 18611 is that baryonic matter at sub¬ 
saturation densities presents a first order liquid-gas phase tran¬ 
sition which is signaled by the fact that a huge part of the 
phase diagram is jumped over if one imposes constant chemi¬ 
cal potentials lf5lll86fl . 

As we have discussed in section iHEl this instability is not 
physical and only comes from the fact that the electron con¬ 
tribution is neglected in the instability analysis. If the electron 
free energy is accounted for, the dependence of the free energy 
density on the baryonic density reads 

f{pB,Pp)=fB{pB,Pp)+fel{pp), 


( 122 ) 



































23 



-8 -7 -6 -5 -4 -3 -2 

10 10 10 10 10 10 10 




1 

10 

10 

10 

10 

10 

10 


10 

10 

10 

10 

10 

10 


10 

10 

10 

10 



10' 8 10 7 10' 6 10' 5 IQ ' 4 10' 3 10' 2 
Pb ( fm ) 



-8 -7 -6 -5 -4 -3 -2 

10 10 10 10 10 10 10 




-7 -6 -5 -4 -3 -2 

10 10 10 10 „ 10 10 
— 0.4 Pb ( fm ) 


0.6 
0.8 
1.0 
1.2 
1.4 
1.6 
1.8 

- 2.0 



FIG. 16. (Color online). Behavior of the constrained free energy 
in j6-equilibrium {jli = 0) for the clusterized phase (full line) and 
the homogeneous phase (dashed line) at different values of the bary- 
onic chemical potential Xb mentioned in the key legend (in MeV). 
Experimental binding energies and predictions of the 10-parameter 
mass model of Duflo-Zuker are used for the nuclear clusters. For the 
unbound nucleon gas the SLY4 effective interaction has been used. 
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FIG. 15. (Color online). NSE mass fractions of unbound nucleons 
and 2 H, 3 H, 3 He, 4 He, A - 4 H and A - 5 He at /3—equilibrium for dif¬ 
ferent densities and temperatures expressed in MeV and listed in the 
key legend. Experimental El and DZ10 (9li] data have been used 
for the binding energies. The SLY4 effective interaction was used for 
the unbound nucleon component. Note that X-axis range is not the 
same in all panels. 

where fn denotes the baryonic part and the Coulomb interac¬ 
tion part between protons and electrons, and f e i = - T I n z r p . 
The relations ( |43| > between density and chemical potential are 
shifted because of the electron contribution 

He -f Mb = Mb + ^Pei \ M3 -f M/ = M3 - ^Peh (123) 

and the curvature of the constrained free energy density is aug¬ 
mented of a positive term as: 

pr = p + ipL ( 124 ) 

opjj "Pb 2 dp e i 

This quenching of the phase transition has as a practical 
consequence that a one-to-one correspondence between den¬ 
sity and chemical potential exists in stellar matter, meaning 


FIG. 17. (Color online). Fragment mass distributions correspond¬ 
ing to /3 -equilibrium, T =2 MeV and different baryonic densities (ex¬ 
pressed in fm~ 3 ) as listed in the key legend. The numbers next to 
peaks specify the charge number of the most abundant nucleus. 


that it is possible to describe all the possible density configura¬ 
tions in a grandcanonical treatment, provided the electron con¬ 
tribution is accounted for. In that case, the ensemble equiva¬ 
lence is recovered and the associated partitions are by con¬ 
struction identical to the ones obtained in a canonical model, 
as we have explicitly demonstrated in sections fill B IIII Cl 

It is however in principle perfectly possible that a resid¬ 
ual convexity persists in the constrained free energy (1 1 24b . In 
that case, a first order phase transition reminiscent of liquid- 
gas would survive in stellar matter. Such a phenomenology 
was recently suggested in ref. l l8l ill and evidenced, for T=0, in 
Fig.6(b). 

To answer to this question, we show in Fig[T6]the compar¬ 
ison between the constrained free energy energy density of 
the clusterized Wigner-Seitz cell, and the one corresponding 
to homogeneous matter, for different values of the baryonic 
chemical potential (mentioned in the key legend), at a repre¬ 
sentative temperature of 2 MeV. We can see that the cluster¬ 
ized phase systematically presents a lower free energy den- 
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sity than the homogeneous system, for all chemical potentials 
up to about j Ub = 14 MeV. For the highest considered chem¬ 
ical potential, 15 MeV, the constrained energy minimum cor¬ 
responds to homogeneous matter. This means that the first or¬ 
der phase transition is restricted to a density domain between 
about Pb = 0.07 and p« = 0.09 fm 1 2 3 4 * . These values obvi¬ 
ously depend on the temperature and on the effective interac¬ 
tion, but still the associated density discontinuity is too small 
to have any observable effects. Moreover, we have to stress 
again that we have disregarded deformation degrees of free¬ 
dom in this model. The inclusion of highly deformed pasta 
clusters would lead to a lowering of the clusterized phase, and 
an extra shrinking of the possible transition domain. 

Figure[l7]shows the detailed matter composition in the high 
density region close to the transition to homogeneous matter. 
Dominance of exotic light nuclei as 7 H, 14 He, 17 Li, 20 Be, 22 Be 
is worthwhile to note meaning that it is very important to ac¬ 
count for light clusters in that domain. It is therefore possible 
that smoother transitions would be observed between the dif¬ 
ferent pasta phases, and between the pasta phase and homo¬ 
geneous matter, if light clusters were accounted for. We leave 
this point to future developments. 

IV. CONCLUSIONS 

In this paper we have presented a unified treatment of the 
stellar matter composition and equation of state in the sub¬ 
saturation regime, which can be applied at any temperature, 
density and proton fraction. 

The basic idea of the model is to consider stellar matter as a 
statistical mixing of independent Wigner-Seitz cells. The in¬ 
dividual composition in terms of bound and unbound particles 
does not minimize the free energy density, but the combina¬ 
tion of different cells does. 

The result is a set of NSE-like equations for the cluster 
abundancies, where both the bulk and the surface part of the 
cluster self-energies are modified by the presence of free nu¬ 
cleon scattering states, and a high energy cut naturally ap¬ 
pears in the cluster internal state partition sum. The model 
dependence of the finite temperature model is thus limited to 
the model dependence of the treatment of the Wigner-Seitz 
cell, which in turn is very well constrained by microscopic 
calculations, with a residual uncertainty limited to the den¬ 
sity dependence of the symmetry energy in the underlying ef¬ 
fective interaction, and the detailed treatment of the isospin 
dependent surface tension. In the present applications, the 
in-medium modifications are treated in the local density ap¬ 


proximation, but it will be extremely interesting to map them 
from a more sophisticated microscopic, and possibly beyond 
mean-field treatment in the next future. 

We have analytically shown that, for a given set of chemical 
potentials, the most probable cell composition coincides with 
the one which is obtained by the standard variational proce¬ 
dure assuming one single representative cluster. This guaran¬ 
tees that the model has the correct zero temperature limit. 

However, the simultaneous presence of many different clus¬ 
ters in each thermodynamic condition modifies the relation 
between density and chemical potential with respect to the sin¬ 
gle nucleus approximation. As a consequence, stellar matter 
predictions of this improved NSE model differ from the SNA 
approximation even at the level of the most probable compo¬ 
sition, and even at temperatures lower than 1 MeV. 

We have specifically shown quantitative applications in /3- 
equilibrium. The dominant configuration is a mixture of clus¬ 
ters of different mass and atomic numbers. This effect is due 
at low density to the non-monotonic behavior of the cluster en¬ 
ergies due to shell and sub-shell closures, and at high density 
to the flatness or multi-minima of the free-energy landscape 
for very neutron rich matter. None of these features can be 
accounted in a SNA approach. In addition to the multi-peaked 
cluster distribution we have seen that very light clusters appear 
with a probability compared to the one of heavier clusters. 
This feature is accounted in the Lattimer-Swesty SNA model 
by including a-particles in the single representative Wigner- 
Seitz cell. We can see that this is physically correct at the 
lowest densities, which at j3-equilibrium correspond to matter 
close to isospin symmetry. Conversely, in very asymmetric 
matter as it can be found at (5 -equilibrium at higher density, 
the most probable light cluster is never a a particle, but rather 
the last bound isotope of H and He. It is therefore clear that 
at finite temperature other light particles than a have to be 
included in the equilibrium. 
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